Chapter 10 The Covariance Approach to Recognition:

How can we identify a person we are looking at with a video camera quickly? This has a lot of applications to maintaining secure access to an industrial site among other things. The technology here is based on what are called eigenfaces , a nice use of a thing called eigenvectors . We can also use these ideas to introduce the idea of searching a huge collection of images (such as NASA has stored at various sites) by the use of what might be called eigenimages .

10.1 The Basic Plan:

The idea here is pretty straightforward. Here are the steps:

10.1.1 Collect Images:

For the MTH SC 450 class during the Fall of 2000, we took 52 pictures of the people in the class. There were 11 students in the class and we added pictures of the instructor Jim, his wife Pauli and daughter Qaitlin. We took 5 pictures per person. We used 4 of these to create what might be called an average image of each person. We tried to make these first 4 pictures sort of uniform. Our approach was to take each picture from the same distance with the class room wall ( a nicely uniform blah color ) as the background. Each of the 4 pictures were supposed to be a bit different; lookng straightforward but with different facial expressions and so forth. The fifth picture was to be something different. Here each person looked to the side or something so that the fifth image was very distinct from the average of the first four.

We got all the pictures on a Kodak PhotoCD which turned out to be jpeg files that were fairly big, 1500 by 1000. So we preprocessed these in various ways to get images we would use in our recognition algorithm.

10.1.2 The Original Pauli Images:

For example, here are the original four images of the domestic engineer and goddess Pauli: these can be seen in Figure 10.1, Figure 10.2, Figure 10.3 and Figure 10.4.

Figure 10.1: Pauli Image 1

Figure 10.2: Pauli Image 2

Figure 10.3: Pauli Image 3

Figure 10.4: Pauli Image 4

Now the images you see are rotated -90 degrees from what was supplied on the CDROM and then resized to 1/8 of the original size. This resulted in a file size of 192 by 128. So how do we do this in Matlab? First, we would read in an image like so:

%read in the jpeg files from disk
>>I1 = imread('image1.jpg');
>>I2 = imread('image2.jpg');
>>I3 = imread('image3.jpg');
>>I4 = imread('image4.jpg');

We would then rotate the images and resize them, but the details of that code are important right now. So let's just assume the images we read in are the right size and orientation for now.

This reads the jpeg files into internal Matlab representations which consist of 192 by 128 three dimensional row vectors of the form [r, g, b]. The r, g and b values are of data type uint8 which is what is known as an unsigned integer of size 8. This means that an r value is represented by a string of 8 binary bits of the form [b₁,b₂,b₃,b₄,b₅,b₆,b₇,b₈] where each bit b_i can be 0 or 1. This bit string is then interpreted as

b₁*2⁷ + b₂*2⁶ + b₃*2⁵ + b₄*2⁴ + b₅*2³ + b₆*2² + b₇*2¹ + b₈

or vice-versa depending on how the computer architecture is set up.

This means that the smallest value is 0 and the largest is 255. The r is called the red value; the g is the green value and the last number of the triple, b is the blue value. The various combinations of these integers are then mapped into colors using a variety of color models. We thus call this the RGB color model for obvious reasons. The triple (0,0,0) is usually black and the triple (255,255,255) is white. For example, the entry at I1(42,67) might be the vector [217, 37, 98].

10.1.3 Computing the Pauli Average Image:

Matlab does not allow us to add data of type unit8 so to compute the average image, we convert the [r,g,b] values in the image matrices to double, manipulate as we desire and then convert back to uint8 to display. The Matlab code might be

% convert to double
>>ID1 = double(I1);
>>ID2 = double(I2);
>>ID3 = double(I3);
>>ID4 = double(I4);

% compute average
>> AV = (ID1+ID2+ID3+ID4)/4.0;

To display this average, we would type:

image( uint8(AV) );

The average of the domestic goddess images can be seen in Figure 10.5.

Figure 10.5: The Pauli Image

10.1.4 The Pauli Test Image:

Now the fifth image will not be used in the computation of an image and we will see later how we use this to test our recognition engine's accuracy.

Figure 10.6: The Pauli Test Image

10.2 Image Processing:

The 13 different people with their 4 images apiece thus give us 52 different images which are converted by the process above to 3 dimensional matrices of double. So there are a lot of numbers for each image!

Call these prepocessed images A_i.

10.2.1 Compute the Average of All Images:

We just compute the average of all the images as follows:

p=1

2 A_p

This image looks like Figure 10.7:

Figure 10.7: The Average Image

10.2.2 Compute Difference Images:

Next, we compute the difference image

D_i

A_i - S

10.2.3 Convert Image Matrices to Vectors:

Now these difference images consist of [r,g,b] triples of double each of which lie between -.75 × 255 and .75 × 255. So it is possible for some of the values in these difference triples to be negative. We will map these raw differences to the interval [0, 255] via

[r,g,b]

([r,g,b] + .75 × 255) ×

10.2.4 Convert The Difference Images to Base 255 Doubles:

We will convert these matrices of [r,g,b] type data to vectors as follows: first convert to a base 255 number:

[r,g,b]

r × 255² + g × 255 + b

Next, scale this to be in [[0,1]:

[r,g,b]

[r,g,b] ×

255³

We call these new matrices of double only M_i.

10.2.5 Convert the Double Matrices to Vectors:

We convert the matrix p by q matrix M to a vector V of size p × q by 1 as follows: for each row i and column j, we set

V_{(i-1)q + j}

M_i,j

Each of the matrices M_i is converted in this fashion into a vector d_i.

These vectors are quite large!

10.2.6 Compute the Covariance Matrix:

Now find the covariance matrix

A A^T

where A is the 192 × 128 by 52 matrix

[d₁ | d₂ | ··· | d₅₂]

This makes the matrix A A^T 192 × 128 by 192 × 128 -- much to large for ease of computation!

10.2.7 Compute Eigenvectors of A^T A:

Theory tells us that the eigenvectors and eigenvalues of the covariance matrix are important for image recognition. However, finding these things for such a huge matrix is way too hard! It turns out the eigenvalues and eigenvectors of the smaller 52 by 52 matrix A^T A are just what we need. There are at most 52 of these and theory says the eigenvalues are nonnegative and the eigenvectors are 90 degrees apart. If we call these eigenvectors

v₁ , ... , v₅₂

then the vectors

A v₁ , ... , A v₅₂

turn out to be eigenvectors for the covariance matrix. These are called eigenfaces .

Now for all of our 52 images, we can do these computations and we can see the first ten resulting eigenfaces in Figure 10.8 to Figure 10.17.

Figure 10.8: EigenImage 1

Figure 10.9: EigenImage 2

Figure 10.10: EigenImage 3

Figure 10.11: EigenImage 4

Figure 10.12: EigenImage 5

Figure 10.13: EigenImage 6

Figure 10.14: EigenImage 7

Figure 10.15: EigenImage 8

Figure 10.16: EigenImage 9

Figure 10.17: EigenImage 10

10.2.8 Find Feature Vectors:

Now project each original difference vector d_i onto the subspace determined by a small number of these eigenfaces. For example, difference image d₃ might look like

d₃

.22*Av₁ + .34*Av₂ -.67*Av₃ + 46.7*Av₄ + Z

where the first four pieces is the portion of d₃ that lies in the plane determined by the first four eigenfaces and Z is the part of d₃ that lies out of this plane. So we could associate d₃ with a small vector called a feature vector:

.22

.34

-.67

46.7

We can compute feature vectors for all of the 52 difference vectors for a given number of eigenvectors. For our problem of 52 images, we have usually chosen 25 eignevectors. Note in effect, we are assigning a small 25 dimensional vector to each large original image.

10.2.9 Compute Class Vectors:

Take the four feature vectors for each person and average them. This average feature vector is called the class vector for that face.

10.2.10 Compute The Test Feature Vector:

Now take a new picture of a person in the original group and convert it to a difference vector like usual.

For example, the test image of the domestic engineer and goddess Pauli is given in Figure 10.6 as previously mentioned. Then find its feature vector and calculate the the distance from this feature vector to all the class vectors. If it is sufficiently close to a class vector, we can say it is an image from the image class for the person the class represents.

For the Pauli test image, the projection of its difference vector to the eigenface subspace can be seen in Figure 10.18.

Figure 10.18: The Projected Pauli Test Image

Now this is just a start at this big subject of facial recognition. We can also use artificial neural technology and other things. This technique is statistically based and you really need some tools like Matlab to pull it off!

Also, how would we choose feature vectors so that we can describe all of the images in the NASA database in such a way that we could use search engines? We would like to be able to query the data base like this: Find us all images of the planet mars that have a spaceship like object in orbit . We have very little clue how to do this but techniques like the one discussed above and many others are being actively considered.

10.3 Implementing The Basic Plan in Matlab:

Now let's begin the task of putting together our implementation of the Covariance approach to Object Recognition using Matlab.

Currently, our implementation requires the following scripts:

Basic Input/ Output Operations
- : ReadDoubleImages.m Reads N by M by 3 (rgb) data from a file
- : WriteDoubleImages.m Writes N by M by 3 rgb data to a file.
Prepocessing and Getting the Images from Disk
- : ProcessFiles.m Reads original images and processes them into the images we want to use to build our recognition engine.
Processing the Images for Recognition Purposes
- : ProcessAllDirectories.m This looks at all the processed image data and computes the auxiliary data structures we need to build the recognition engine. This code uses the following functions:
  - : GetDirectoryAverage.m Find the average of the image data in each directory.
  - : ComputeAverage.m Get the average of all the image data.
  - : GetDiffs.m Find the difference between the image data and the average of the image data.
  - : rgb2double.m Convert the differenced image data from the rgb format to a single double at each pixel position.
  - : matrix2vector.m Conver the matrix above into a vector.
Finding the EigenImages
- : GetEigenImages.m This finds the eigenimages for our data.
The Class and Feature Vectors
- : GetClassVectors.m Each image class has featrue vectors. We average these to get the Class feature vector.
- : GetFeatureVectors.m Each of the images will have a feature vector. This functions computes all of them.
The Full Build Process:
- : BuildCovariance.m Build the Covariance Matrix using our scripts.
- : BuildRecognition.m Build the Recognition Engine using the Covariance matrix.
Testing the Recognition Engine:
- : TestImage.m We read in a test image and compute its test feature. We use the functions:
  - : GetTestFeature.m This computes the test feature vector for the test image.
  - : getTestDistance.m This finds the distance of the test feature vector from each class feature vector.
Conversion of Vector Data Back to Image Data
- : vector2matrix.m This converts the vector data back to an image.

10.4 Using the Functions In a Matlab Session:

A typical Matlab session is as follows:

% let's assume you already have the image files in 
% the right size stored on disk
%

%
%  set up paths
%
>> path(path,'/local2/petersj/eigenfaces/Images');

If not, we would have called the function PrcessFiles with code:

ProcessFiles.m

 1 function ProcessFiles(name)
 2 %
 3 % process files in given directory name
 4 %
 5 
 6 %
 7 %
 8 %
 9 cd(name);
10 
11 %
12 % read in full images
13 %
14 I1 = imread('image1.jpg');
15 I2 = imread('image2.jpg');
16 I3 = imread('image3.jpg');
17 I4 = imread('image4.jpg');
18 %
19 % rotate images -90 degrees
20 %
21 IR1 = imrotate(I1,-90);
22 IR2 = imrotate(I2,-90);
23 IR3 = imrotate(I3,-90);
24 IR4 = imrotate(I4,-90);
25 %
26 % resize images to 1/8 full size
27 %
28 IR1125 = imresize(IR1,.125);
29 IR2125 = imresize(IR2,.125);
30 IR3125 = imresize(IR3,.125);
31 IR4125 = imresize(IR4,.125);
32 %
33 % write out resized images to files
34 %
35 imwrite(IR1125,'image1125.jpg');
36 imwrite(IR2125,'image2125.jpg');
37 imwrite(IR3125,'image3125.jpg');
38 imwrite(IR4125,'image4125.jpg');
39 
40 cd('..');

syntax highlighted by Code2HTML, v. 0.8.8b

We then need to process things in the directories that store the files. We will use this function: function Av = ProcessAllDirectories which has code:

ProcessAllDirectories.m

 1 function Av = ProcessAllDirectories
 2 
 3 p = 192;
 4 q = 128;
 5 Av = zeros(p,q,3,13);
 6 
 7 B = GetDirectoryAverage('bayless');
 8 Av(1:p,1:q,1:3,1) = B;
 9 
10 B = GetDirectoryAverage('cole');
11 Av(1:p,1:q,1:3,2) = B;
12 
13 B = GetDirectoryAverage('digiuseppe');
14 Av(1:p,1:q,1:3,3) = B;
15 
16 B = GetDirectoryAverage('elmore');
17 Av(1:p,1:q,1:3,4) = B;
18 
19 B = GetDirectoryAverage('fawcett');
20 Av(1:p,1:q,1:3,5) = B; 
21 
22 B = GetDirectoryAverage('harris');
23 Av(1:p,1:q,1:3,6) = B;
24 
25 B = GetDirectoryAverage('jim');
26 Av(1:p,1:q,1:3,7) = B;
27 
28 B = GetDirectoryAverage('kashema');
29 Av(1:p,1:q,1:3,8) = B;
30 
31 B = GetDirectoryAverage('mcleod');
32 Av(1:p,1:q,1:3,9) = B;
33 
34 B = GetDirectoryAverage('pauli');
35 Av(1:p,1:q,1:3,10) = B;
36 
37 B = GetDirectoryAverage('qaitlin');
38 Av(1:p,1:q,1:3,11) = B;
39 
40 B = GetDirectoryAverage('terry');
41 Av(1:p,1:q,1:3,12) = B;
42 
43 B = GetDirectoryAverage('wilkins');
44 Av(1:p,1:q,1:3,13) = B;

syntax highlighted by Code2HTML, v. 0.8.8b

This uses the function function Av = GetDirectoryAverage(name) to do the work which has source:

GetDirectoryAverage.m

 1 function Av = GetDirectoryAverage(name)
 2 %
 3 % get average of images in directory
 4 %
 5 
 6 %
 7 % go to directory
 8 %
 9 cd(name);
10 
11 %
12 %  read in rotated 1/8 scale images
13 %
14 IR1125 = imread('image1125.jpg');
15 IR2125 = imread('image2125.jpg');
16 IR3125 = imread('image3125.jpg');
17 IR4125 = imread('image4125.jpg');
18 
19 %
20 % convert to doubles
21 %
22 ID1 = double(IR1125);
23 ID2 = double(IR2125);
24 ID3 = double(IR3125);
25 ID4 = double(IR4125);
26 
27 %
28 % write double images to files
29 %
30 fd = fopen('dimage1','w');
31 fwrite(fd,ID1,'double');
32 fclose(fd);
33 
34 fd = fopen('dimage2','w');
35 fwrite(fd,ID2,'double');
36 fclose(fd);
37 
38 fd = fopen('dimage3','w');
39 fwrite(fd,ID3,'double');
40 fclose(fd);
41 
42 fd = fopen('dimage4','w');
43 fwrite(fd,ID4,'double');
44 fclose(fd);
45 
46 %
47 % Find average
48 %
49 Av = (ID1+ID2+ID3+ID4)/4.0;
50 
51 cd('..');

syntax highlighted by Code2HTML, v. 0.8.8b

This function will compute the double versions of the images and save them to disk as dimage<k> for k 1:4 and compute the average and return it to the matlab session. So to do all this we just call the ProcessAllDirectories function:

%
% Process all directories, save double versions to
% disk and compute averages and store in 
% Av which is 192 x 128 x 3 x 13
%

>> Av = ProcessAllDirectories

And get the average by calling ComputeAverage which finds the average for all the images. This function has the source:

ComputeAverage.m

 1 function Average = ComputeAverage(Av)
 2 
 3 %
 4 % Av is 192 x 128 x 3 x 13
 5 %
 6 [p,q,r,s] = size(Av);
 7 Average = zeros(p,q,r);
 8 for count = 1:s
 9   Average = Average + Av(1:p,1:q,1:3,count);
10 end
11 Average = Average/s;

syntax highlighted by Code2HTML, v. 0.8.8b

The matlab call is

>> average = ComputeAverage;
\end{verbatgim}

If we want to, we can then display our average by typing
 
\begin{verbatim}
%
% Display average
%
>> image( uint8(average));

For later use, we can write this to disk with these lines:

%
% write average to disk
%
fd = fopen('average','w');
fwrite(fd,average,'double');
fclose(fd);

Next, we need to compute the difference images and convert them to vectors: we do this with a call to the function function [diff,v]= GetDiffs(name,S,Average) which has code:

GetDiffs.m

 1 function [diff,v]= GetDiffs(name,S,Average)
 2 %
 3 % get image-average of images in directory
 4 %
 5 
 6 %
 7 % go to directory
 8 %
 9 cd(name);
10 
11 %
12 %  read in rotated 1/8 scale images
13 %
14 IR1125 = imread('image1125.jpg');
15 IR2125 = imread('image2125.jpg');
16 IR3125 = imread('image3125.jpg');
17 IR4125 = imread('image4125.jpg');
18 
19 %
20 % convert to doubles
21 %
22 ID1 = double(IR1125);
23 ID2 = double(IR2125);
24 ID3 = double(IR3125);
25 ID4 = double(IR4125);
26 
27 %
28 % Find differences
29 %
30 [p,q,s] = size(ID1);
31 r = p*q;
32 diff = zeros(p,q,s,4);
33 diff(1:p,1:q,1:3,1) = ID1-Average;
34 diff(1:p,1:q,1:3,2) = ID2-Average;
35 diff(1:p,1:q,1:3,3) = ID3-Average;
36 diff(1:p,1:q,1:3,4) = ID4-Average;
37 
38 %
39 % convert rgb triples to doubles
40 % now the diff's are in the range (-.75*255, .75*255)
41 % the functoion rgb2double expects to see an
42 % input with entries in the range [0,1]
43 % so we process as 
44 % (diff(1:p,1:q,1:3,*)+ .75*255)*(2.0/3.0)
45 %
46 m1 = rgb2double((diff(1:p,1:q,1:3,1)+ .75*255)*(2.0/3.0));
47 m2 = rgb2double((diff(1:p,1:q,1:3,2)+ .75*255)*(2.0/3.0));
48 m3 = rgb2double((diff(1:p,1:q,1:3,3)+ .75*255)*(2.0/3.0));
49 m4 = rgb2double((diff(1:p,1:q,1:3,4)+ .75*255)*(2.0/3.0));
50 
51 %
52 % convert to vectors
53 %
54 [p,q] = size(m1);
55 r = p*q;
56 v = zeros(r,4);
57 v(1:r,1) = transpose(matrix2vector(m1));
58 v(1:r,2) = transpose(matrix2vector(m2));
59 v(1:r,3) = transpose(matrix2vector(m3));
60 v(1:r,4) = transpose(matrix2vector(m4));
61 
62 cd('..');

syntax highlighted by Code2HTML, v. 0.8.8b

The Matlab lines to do this are do directory by directory:

%
% so get all diffs
%
>> S = size(average);
>> [diff1,v1]= GetDiffs('bayless'   ,S,average   );
>> [diff2,v2]= GetDiffs('cole'      ,S,average   );
>> [diff3,v3]= GetDiffs('digiuseppe',S,average   );
>> [diff4,v4]= GetDiffs('elmore'    ,S,average   );
>> [diff5,v5]= GetDiffs('fawcett'   ,S,average   );
>> [diff6,v6]= GetDiffs('harris'    ,S,average   );
>> [diff7,v7]= GetDiffs('jim'       ,S,average   );
>> [diff8,v8]= GetDiffs('kashema'   ,S,average   );
>> [diff9,v9]= GetDiffs('mcleod'    ,S,average   );
>> [diff10,v10]= GetDiffs('pauli'     ,S,average   );
>> [diff11,v11]= GetDiffs('qaitlin'   ,S,average   );
>> [diff12,v12]= GetDiffs('terry'     ,S,average   );
>> [diff13,v13]= GetDiffs('wilkins'   ,S,average   );

Next, we set up the matrix of difference vectors T:

%
% Set Up Difference Matrix
%
>> [p,q] = size(v1);
>> T = zeros(p,52);
>> T(1:p,1:4) = v1;
>> T(1:p,5:8) = v2;
>> T(1:p,9:12) = v3; 
>> T(1:p,13:16) = v4; 
>> T(1:p,17:20) = v5;
>> T(1:p,21:24) = v6; 
>> T(1:p,25:28) = v7; 
>> T(1:p,29:32) = v8; 
>> T(1:p,33:36) = v9; 
>> T(1:p,37:40) = v10; 
>> T(1:p,41:44) = v11; 
>> T(1:p,45:48) = v12; 
>> T(1:p,49:52) = v13;

Now we can compute the covariance matrix itself:

%
% Find Covariance
%
>> TT = transpose(T);
>> W = (TT*T)/52.0;

We can then write it to disk:

%
% write covariance to disk
%
fd = fopen('covariance','w');
fwrite(fd,W,'double');
fclose(fd);

Note if we wanted to read it back, we would type

fd = fopen('covariance','r');
S = size(W);
W = fread(fd,S,'double');

Then we find the first 20 eigenvalues and associated eigenvectors:

>> [EigenVs,EigenVals] = eigs(W,20);
>> EigenVals

EigenVals =

   1.0e+03 *

  Columns 1 through 7 

    6.2256         0         0         0         0         0         0
         0    0.1008         0         0         0         0         0
         0         0    0.0350         0         0         0         0
         0         0         0    0.0226         0         0         0
         0         0         0         0    0.0173         0         0
         0         0         0         0         0    0.0109         0
         0         0         0         0         0         0    0.0102
         0         0         0         0         0         0         0
         0         0         0         0         0         0         0
         0         0         0         0         0         0         0
         0         0         0         0         0         0         0
         0         0         0         0         0         0         0
         0         0         0         0         0         0         0
         0         0         0         0         0         0         0
         0         0         0         0         0         0         0
         0         0         0         0         0         0         0
         0         0         0         0         0         0         0
         0         0         0         0         0         0         0
         0         0         0         0         0         0         0
         0         0         0         0         0         0         0

  Columns 8 through 14 

         0         0         0         0         0         0         0
         0         0         0         0         0         0         0
         0         0         0         0         0         0         0
         0         0         0         0         0         0         0
         0         0         0         0         0         0         0
         0         0         0         0         0         0         0
         0         0         0         0         0         0         0
    0.0086         0         0         0         0         0         0
         0    0.0073         0         0         0         0         0
         0         0    0.0065         0         0         0         0
         0         0         0    0.0055         0         0         0
         0         0         0         0    0.0052         0         0
         0         0         0         0         0    0.0048         0
         0         0         0         0         0         0    0.0044
         0         0         0         0         0         0         0
         0         0         0         0         0         0         0
         0         0         0         0         0         0         0
         0         0         0         0         0         0         0
         0         0         0         0         0         0         0
         0         0         0         0         0         0         0

  Columns 15 through 20 

         0         0         0         0         0         0
         0         0         0         0         0         0
         0         0         0         0         0         0
         0         0         0         0         0         0
         0         0         0         0         0         0
         0         0         0         0         0         0
         0         0         0         0         0         0
         0         0         0         0         0         0
         0         0         0         0         0         0
         0         0         0         0         0         0
         0         0         0         0         0         0
         0         0         0         0         0         0
         0         0         0         0         0         0
         0         0         0         0         0         0
    0.0040         0         0         0         0         0
         0    0.0037         0         0         0         0
         0         0    0.0036         0         0         0
         0         0         0    0.0034         0         0
         0         0         0         0    0.0032         0
         0         0         0         0         0    0.0029

Then we can compute the EigenImages with a call to function EI = GetEigenImages(n,T,EigenVs) with code:

GetEigenImages.m

 1 function EI = GetEigenImages(n,T,EigenVs)
 2 
 3 [p,q] = size(T);
 4 [r,s] = size(EigenVs);
 5 EI = zeros(p,n);
 6 for count = 1:n
 7   A = T*EigenVs(1:r,count);
 8   length = norm(A);
 9   EI(1:p,count) = A/length;
10 end

syntax highlighted by Code2HTML, v. 0.8.8b

We type the following Matlab lines:

>> EI = GetEigenImages(10,T,EigenVs);

Next, we get the feature vectors by calling function F = GetFeatureVectors(n,v,EI) which has code:

GetFeatureVectors.m

1 function F = GetFeatureVectors(n,v,EI)
2 
3 F = zeros(n,4);
4 [p,q] = size(EI);
5 for count = 1:4
6   for row = 1:n
7       F(row,count) = transpose(v(1:p,count))*EI(1:p,row);
8     end
9 end

syntax highlighted by Code2HTML, v. 0.8.8b

In Matalb we type:

>> F = zeros(20,52);
>> F(1:20,1:4) = GetFeatureVectors(20,v1,EI);
>> F(1:20,5:8) = GetFeatureVectors(20,v2,EI); 
>> F(1:20,9:12) = GetFeatureVectors(20,v3,EI);
>> F(1:20,13:16) = GetFeatureVectors(20,v4,EI); 
>> F(1:20,17:20) = GetFeatureVectors(20,v5,EI);
>> F(1:20,21:24) = GetFeatureVectors(20,v6,EI);
>> F(1:20,25:28) = GetFeatureVectors(20,v7,EI); 
>> F(1:20,29:32) = GetFeatureVectors(20,v8,EI);
>> F(1:20,33:36) = GetFeatureVectors(20,v9,EI); 
>> F(1:20,37:40) = GetFeatureVectors(20,v10,EI);
>> F(1:20,41:44) = GetFeatureVectors(20,v11,EI); 
>> F(1:20,45:48) = GetFeatureVectors(20,v12,EI);  
>> F(1:20,49:52) = GetFeatureVectors(20,v13,EI);

So for example the feature vector for the first image in the first directory is

>> F(1:20,1)

ans =

  -85.8201
   -9.7383
   -1.3781
    4.3543
    0.5797
    3.0650
    1.0536
    2.2167
   -1.0532
   -0.1967
   -1.0027
    0.1949
    0.0793
    0.3115
    0.2228
   -1.9004
   -1.0676
    0.2727
    1.2976
   -0.6332

In this work we also use some other utility functions: matrix2vector, rgb2double and vector2matrix with codes

matrix2vector.m

 1 function V = matrix2vector(A)
 2 %
 3 % A is N x M matrix of doubles
 4 % we will convert as follows:
 5 % A(p,q) ==> V((p-1)*M + q
 6 % So first row of A goes to V(1) ...V(M)
 7 %    second row of A goes to V((M+1) to V(2M)
 8 %    last row of A goes to V( (N-1)*M ) to V(N*M)
 9 %
10 [N,M] = size(A);
11 size(A);
12 for p=1:N
13   for q=1:M
14     k = (p-1)*M+q;
15     V(k) = A(p,q);
16   end
17 end

syntax highlighted by Code2HTML, v. 0.8.8b rgb2double.m

 1 function B = rgb2double(A)
 2 %
 3 % A is M x N x 3 data
 4 % Each entry in A is a rgb triple, [r, g, b]
 5 % r, g and b are in [0,255]
 6 % The result is a M x N matrix B
 7 % whose entries are doubles given by
 8 % r* 256^2 + g*256 + b
 9 % which can be done more efficiently as
10 % 256*(256*r + g) + b
11 %
12 [M,N,P] = size(A);
13 R = 256*256*256;
14 B = zeros(M,N);
15 size(B);
16 for j=1:N
17   B(1:M,j) =  ( 256*( A(1:M,j,1)*256+A(1:M,j,2) ) + A(1:M,j,3))/R;
18 end

syntax highlighted by Code2HTML, v. 0.8.8b vector2matrix.m

 1 function A = vector2matrix(V,N,M)
 2 %
 3 % V is a N*M vector which will be
 4 % converted back to a N x M matrix of 
 5 % [r g b] values
 6 % we will convert as follows:
 7 % for index p in V, 
 8 %   V(p) =  (r*256^2 + g*256 + b)/256^3
 9 %   scale by 256^3 giving
10 %   x = r*256^2 + g*256 + b
11 %   compute 
12 %   y = x/256 = r*256 + g + fraction b/256
13 %   convert to int
14 %   yint = double(int16(y)) = r*256+g
15 %   compute 
16 %   b = x - 256*yint
17 %   compute
18 %   z = y/256 = r + g/256 + b/256^2
19 %   convert to int 
20 %   r = double(int16(z))
21 %   compute
22 %   g = yint - 256*r 
23 %
24 P = N*M;
25 A = zeros(N,M,3);
26 for n=1:N
27   for m = 1:M
28     p = (n-1)*M+m;
29     % V(p) is a double
30     x = V(p)*256*256*256;
31     y = x/256;
32     yint = double(int32(y));
33     b = x - 256*yint;
34     z = y/256;
35     r = double(int32(z));
36     g = yint - 256*r;
37     A(n,m,1) = r;
38     A(n,m,2) = g;
39     A(n,m,3) = b;
40     %check = (r*256*256+g*256+b)/(256*256*256)
41     %check2 = V(p) - check
42   end   
43 end

syntax highlighted by Code2HTML, v. 0.8.8b

10.5 Using the Build Functions:

%
%  A new set of scripts
%

%
%  set up paths
%
>> path(path,'/local2/petersj/eigenfaces/Images');

%
% Let's semi automate our building of the 
% covariance matrix
%
% we use the function BuildCovariance()
% So to build the covariance matrix, we just call this
% function.  We have put in some printed messages so we
% don't get bored looking at the prompt!!
% 
>> BuildCovariance(192,128);
process all directories
Getting Directory Averages in directory
bayless
Getting Directory Averages in directory
cole
Getting Directory Averages in directory
digiuseppe
Getting Directory Averages in directory
elmore
Getting Directory Averages in directory
fawcett
Getting Directory Averages in directory
harris
Getting Directory Averages in directory
jim
Getting Directory Averages in directory
kashema
Getting Directory Averages in directory
mcleod
Getting Directory Averages in directory
pauli
Getting Directory Averages in directory
qaitlin
Getting Directory Averages in directory
terry
Getting Directory Averages in directory
wilkins
get full average
write average to disk
Get all difference images
Getting Difference Images in directory
bayless
    read in double images from disk
    compute double difference images
    convert double difference images to matrices
    compute matrices to vectors
     write Vectors to disk
Getting Difference Images in directory
cole
    read in double images from disk
    compute double difference images
    convert double difference images to matrices
    compute matrices to vectors
     write Vectors to disk
Getting Difference Images in directory
digiuseppe
    read in double images from disk
    compute double difference images
    convert double difference images to matrices
    compute matrices to vectors
     write Vectors to disk
Getting Difference Images in directory
elmore
    read in double images from disk
    compute double difference images
    convert double difference images to matrices
    compute matrices to vectors
     write Vectors to disk
Getting Difference Images in directory
fawcett
    read in double images from disk
    compute double difference images
    convert double difference images to matrices
    compute matrices to vectors
     write Vectors to disk
Getting Difference Images in directory
harris
    read in double images from disk
    compute double difference images
    convert double difference images to matrices
    compute matrices to vectors
     write Vectors to disk
Getting Difference Images in directory
jim
    read in double images from disk
    compute double difference images
    convert double difference images to matrices
    compute matrices to vectors
     write Vectors to disk
Getting Difference Images in directory
kashema
    read in double images from disk
    compute double difference images
    convert double difference images to matrices
    compute matrices to vectors
     write Vectors to disk
Getting Difference Images in directory
mcleod
    read in double images from disk
    compute double difference images
    convert double difference images to matrices
    compute matrices to vectors
     write Vectors to disk
Getting Difference Images in directory
pauli
    read in double images from disk
    compute double difference images
    convert double difference images to matrices
    compute matrices to vectors
     write Vectors to disk
Getting Difference Images in directory
qaitlin
    read in double images from disk
    compute double difference images
    convert double difference images to matrices
    compute matrices to vectors
     write Vectors to disk
Getting Difference Images in directory
terry
    read in double images from disk
    compute double difference images
    convert double difference images to matrices
    compute matrices to vectors
     write Vectors to disk
Getting Difference Images in directory
wilkins
    read in double images from disk
    compute double difference images
    convert double difference images to matrices
    compute matrices to vectors
     write Vectors to disk
Set up T matrix: holds all difference vectors
write T to disk
Find Covariance
Write covariance to disk

%
%  Next we build the recognition engine
%
>> BuildRecognition(52,192,128,25);
   read in the covariance matrix
   read in the difference vectors

Get

N =

    25

eigenvalues and eigenvectors

eigs =

   1.0e+03 *

    6.2256
    0.1008
    0.0350
    0.0226
    0.0173
    0.0109
    0.0102
    0.0086
    0.0073
    0.0065
    0.0055
    0.0052
    0.0048
    0.0044
    0.0040
    0.0037
    0.0036
    0.0034
    0.0032
    0.0029


stopcrit =

   1.9147e-15

==========================
Get EigenImages
write eigenImages to disk
Get Feature Vectors
write Feature Vectors to disk
Get Class Vectors
write Class Vectors to disk

%
% Now read in average and EigenImages
%

>> fd = fopen('EigenImages');
>> EI = fread(fd,[192*128 25],'double');
>> fclose(fd);
>> size(EI)

ans =

       24576          25

>> fd = fopen('average');
>> fclose(fd);                   
>> average = ReadDoubleImages(192,128,'average');
>> size(average)

ans =

   192   128     3

%
% read in ClassVectors
%
>> fd = fopen('ClassVectors');
>> ClassVectors = fread(fd,[25 13],'double');
>> size(ClassVectors)

ans =

    25    13

%
% Now see if we can classify the test image in  say 'jim'
%
>> dist = TestImage(192,128,'jim','test.jpg',average,EI,ClassVectors);
get a test image from directory
jim
 with name
test.jpg
   rotate image -90 degress
   resize to 1/8 scale
   write new image to disk
  convert image to double
  compute difference image
convert rgb to double
  compute vector
  get test feature vector
   find distance to class vectors
>> dist

dist =

   45.0095
   45.0325
   32.1367
   23.9600
   34.6160
   44.4723
    7.7621
   30.8399
   27.8847
   39.7993
   36.1505
   42.4495
   41.5786

%
% The smallest distance is in component 7.
% Now the order that the directories are
% processed is
%

1.  'bayless'   
2.  'cole'      
3.  'digiuseppe'
4.  'elmore'    
5.  'fawcett'   
6.  'harris'    
7.  'jim'       
8.  'kashema'   
9.  'mcleod'    
10. 'pauli'    
11. 'qaitlin'  
12. 'terry'    
13. 'wilkins' 

%
%  So we classified the test 'jim' image as the 'jim'
%  class vector
%

% 
% Let's try 'kashema'
%

>> dist = TestImage(192,128,'kashema','test.jpg',average,EI,ClassVectors);
get a test image from directory
kashema
 with name
test.jpg
   rotate image -90 degrees
   resize to 1/8 scale
   write new image to disk
   convert image to double
   compute difference image
   convert rgb to double
   compute vector
   get test feature vector
   find distance to class vectors
>> dist

dist =

   25.4525
   25.0242
   16.2876
   21.1934
   17.1000
   24.5143
   27.3750
   12.1808
   15.2702
   20.6883
   20.5079
   20.6519
   20.1377

%
% the smallest is component 8 which is the 'kashema' data!
% Quite Cool!
%

Now here are our functions we have been using:
========================================================================
function BuildCovariance(rows,cols)

%
%  Process all directories
%
disp('process all directories');
Av = ProcessAllDirectories;

%
% Get complete average
%
disp('get full average');
average = ComputeAverage(Av);

%
% Write average to disk
%
disp('write average to disk');
WriteDoubleImages(rows,cols,average,'average');

%
%  Compute difference images
%  Convert to double images in [0 255] range
%  Convert to double 2D matrices in (r*256^2+g*256+b)/256^3 range
%  Convert to 1D vector
%
disp('Get all difference images');
S = size(average);
[diff1,v1]= GetDiffs('bayless'   ,S,average   );
[diff2,v2]= GetDiffs('cole'  ,S,average   );
[diff3,v3]= GetDiffs('digiuseppe',S,average   );
[diff4,v4]= GetDiffs('elmore'  ,S,average   );
[diff5,v5]= GetDiffs('fawcett'   ,S,average   );
[diff6,v6]= GetDiffs('harris'  ,S,average   );
[diff7,v7]= GetDiffs('jim'  ,S,average   );
[diff8,v8]= GetDiffs('kashema'   ,S,average   );
[diff9,v9]= GetDiffs('mcleod'  ,S,average   );
[diff10,v10]= GetDiffs('pauli'     ,S,average );
[diff11,v11]= GetDiffs('qaitlin'   ,S,average );
[diff12,v12]= GetDiffs('terry'     ,S,average );
[diff13,v13]= GetDiffs('wilkins'   ,S,average );

%
% Set Up T Matrix
%
disp('Set up T matrix: holds all difference vectors');
[p q] = size(v1);
T = zeros(p,52);
T(1:p,1:4) = v1;
T(1:p,5:8) = v2;
T(1:p,9:12) = v3; 
T(1:p,13:16) = v4; 
T(1:p,17:20) = v5;
T(1:p,21:24) = v6; 
T(1:p,25:28) = v7; 
T(1:p,29:32) = v8; 
T(1:p,33:36) = v9; 
T(1:p,37:40) = v10; 
T(1:p,41:44) = v11; 
T(1:p,45:48) = v12; 
T(1:p,49:52) = v13; 

%
% write T to disk
%
disp('write T to disk');
fd = fopen('DiffVectors','w');
fwrite(fd,T,'double');
fclose(fd);

%
% Find Covariance
%
disp('Find Covariance');
TT = transpose(T);
W = (TT*T)/52.0;

%
% write covariance to disk
%
disp('Write covariance to disk');
fd = fopen('covariance','w');
fwrite(fd,W,'double');
fclose(fd);

===================================================================
function BuildRecognition(datasize,rows,cols,N)

%
% read in covariance
%
disp('   read in the covariance matrix');
fd = fopen('covariance','r');
covariance = fread(fd,[datasize datasize],'double');
fclose(fd);

%
% read in T
%
disp('   read in the difference vectors');
fd = fopen('DiffVectors','r');
p = rows*cols;
T = fread(fd,[p datasize],'double');
fclose(fd);

%
% Get the desired number of eigenvectors and eigenvalues
%
disp('Get');
N 
disp('eigenvalues and eigenvectors');
[EigenVs,EigenVals] = eigs(covariance,N);

%
% Get EigenImages
%
disp('Get EigenImages');
EI = GetEigenImages(N,T,EigenVs);

%
% Write EigenImages to Disk
%
disp('write eigenImages to disk');
fd = fopen('EigenImages','w');
fwrite(fd,EI,'double');
fclose(fd);

%
% Get Feature Vectors
%
[p q] = size(T);
disp('Get Feature Vectors');
F = zeros(N,datasize);
F(1:N,1:4)   = GetFeatureVectors(N,T(1:p,1:4),EI);
F(1:N,5:8)   = GetFeatureVectors(N,T(1:p,5:8),EI); 
F(1:N,9:12)  = GetFeatureVectors(N,T(1:p,9:12),EI);
F(1:N,13:16) = GetFeatureVectors(N,T(1:p,13:16),EI); 
F(1:N,17:20) = GetFeatureVectors(N,T(1:p,17:20),EI);
F(1:N,21:24) = GetFeatureVectors(N,T(1:p,21:24),EI);
F(1:N,25:28) = GetFeatureVectors(N,T(1:p,25:28),EI); 
F(1:N,29:32) = GetFeatureVectors(N,T(1:p,29:32),EI);
F(1:N,33:36) = GetFeatureVectors(N,T(1:p,33:36),EI); 
F(1:N,37:40) = GetFeatureVectors(N,T(1:p,37:40),EI);
F(1:N,41:44) = GetFeatureVectors(N,T(1:p,41:44),EI);
F(1:N,45:48) = GetFeatureVectors(N,T(1:p,45:48),EI);
F(1:N,49:52) = GetFeatureVectors(N,T(1:p,49:52),EI);

%
% Write Feature Vectors to Disk
%
disp('write Feature Vectors to disk');
fd = fopen('FeatureVectors','w');
fwrite(fd,F,'double');
fclose(fd);

%
% Get Class Vectors
%
disp('Get Class Vectors');
ClassVectors = GetClassVectors(F);

%
% Write Class Vectors to Disk
%
disp('write Class Vectors to disk');
fd = fopen('ClassVectors','w');
fwrite(fd,ClassVectors,'double');
fclose(fd);

=======================================================================
function dist = TestImage(rows,cols,name,filename,average,EI,ClassVectors)

cd(name);
disp('get a test image from directory');
disp(name);
disp(' with name');
disp(filename);
test = imread(filename);

disp('   rotate image -90 degrees');
testrotated = imrotate(test,-90);

disp('   resize to 1/8 scale');
testR125 = imresize(testrotated,.125);

disp('   write new image to disk');
imwrite(testR125,'test125.jpg');

clear test;
clear testrotated;

disp('   convert image to double');
testD = double(testR125);

disp('   compute difference image');
testdiff = ((testD - average)+.75*255)*(2/3);

disp('   convert rgb to double');
testRGB2Double = rgb2double(testdiff);

disp('   compute vector');
testVector = matrix2vector(testRGB2Double);

%
% get test feature vector
%
disp('   get test feature vector');
TestFeature = GetTestFeature(testVector,EI);

 
%
% find distance from our test feature vector 
% to all ClassVectors
%
disp('   find distance to class vectors');
dist = getTestDistance(TestFeature,ClassVectors);

cd('..'); 

===================================================================
function WriteDoubleImages(rows,cols,Data,name)

%
% Write image data that is rows x cols x [rgb] to disk: 
%

%
% Set up temporary matrices to hold R, G, B data
%
R = zeros(rows,cols);
G = zeros(rows,cols);
B = zeros(rows,cols);

R = Data(1:rows,1:cols,1);
G = Data(1:rows,1:cols,2);
B = Data(1:rows,1:cols,3);

fid = fopen(name,'a');
fwrite(fid,R,'double');
fwrite(fid,G,'double');
fwrite(fid,B,'double');

fclose(fid);

================================================================
function A = ReadDoubleImages(rows,cols,name)

%
% Read image data that is in three separate
% sets of data 
% rows x cols for red
% rows x cols for green
% rows x cols for blue
% from disk
%

%
% Set up matrix to hold double image data
%
A = zeros(rows,cols,3);

fid = fopen(name,'r');
A(1:rows,1:cols,1) = fread(fid,[rows cols],'double');
A(1:rows,1:cols,2) = fread(fid,[rows cols],'double');
A(1:rows,1:cols,3) = fread(fid,[rows cols],'double');

fclose(fid);

================================================================

Chapter 10 The Covariance Approach to Recognition:

10.1 The Basic Plan:

10.1.1 Collect Images:

10.1.2 The Original Pauli Images:

10.1.3 Computing the Pauli Average Image:

10.1.4 The Pauli Test Image:

10.2 Image Processing:

10.2.1 Compute the Average of All Images:

10.2.2 Compute Difference Images:

10.2.3 Convert Image Matrices to Vectors:

10.2.4 Convert The Difference Images to Base 255 Doubles:

10.2.5 Convert the Double Matrices to Vectors:

10.2.6 Compute the Covariance Matrix:

10.2.7 Compute Eigenvectors of AT A:

10.2.8 Find Feature Vectors:

10.2.9 Compute Class Vectors:

10.2.10 Compute The Test Feature Vector:

10.3 Implementing The Basic Plan in Matlab:

10.4 Using the Functions In a Matlab Session:

10.5 Using the Build Functions:

10.2.7 Compute Eigenvectors of A^T A: