Your name: Jerry Liu
Your UID: 404474229
Please upload only this notebook to CCLE by the deadline.
Policy for late submission of solutions: We will use Paul Eggert's Late Policy: $N$ days late $\Leftrightarrow$ $2^N$ points deducted} The number of days late is $N=0$ for the first 24 hrs, $N=1$ for the next 24 hrs, etc., and if you submit an assignment $H$ hours late, $2^{\lfloor \,H\,/\,24\, \rfloor}$ points are deducted.
Problem 1: Eigenfaces¶
Chapter 11 of the Course Reader is on Eigenfaces. For this assignment we have included the face files for this chapter in the directory old_faces. It includes some Matlab scripts and a database of 177 face images, each a grayscale .bmp bitmap file of size $64 \times 64$ pixels. The face images have been pre-processed so that the background and hair are removed and the faces have similar lighting conditions.
The Course Reader explains how to reshape each face image into a $1 \times 64^2 \, = \, 1 \times 4096$ row vector, and collect them into a matrix. The principal components of the matrix then define the main dimensions of variance in the faces. The program more_efficient_eigenfaces.m shows how to do this. These principal components are called eigenfaces.
1a: The Average Face¶
Modify the program more_efficient_eigenfaces.m (available in this directory) to use the new_faces images instead of the old_faces images. Also, modify it to use the Matlab function imresize to downsample each of the new faces by a factor of 3, so it is $64 \times 54$ (instead of $193 \times 162$). Then: pad both sides of the image with zeros(64,5) so the result is a $64 \times 64$ image.
Then: create a function that takes as input a string array of filenames of face images, an integer $k$, and an integer sample size $s$ --- and yields the average face and the first $k$ singular values and eigenfaces as output values for a sample of size $s$.
Apply your function to both the new_faces/faces and the new_faces/smiling_faces directories, and plot the absolute value of the difference between your average face and (your downsampled version of) the average face provided in the directory.
(The imagesc function can display images with automatic rescaling of numeric values.)
row = 64;
col = 64;
%----------------------------------------------------------------------------------
% Useful functions for converting images <--> vectors
%----------------------------------------------------------------------------------
%
% These bitmap images are 64x64 matrices of uint8 (unsigned 8-bit) values;
% we convert these to double values in order to permit arithmetic.
%
% We also convert 64 x 64 image matrices to 64^2 x 1 vectors (using reshape())
% so that we can derive eigenfaces using Matlab's SVD/matrix functions.
%
%----------------------------------------------------------------------------------
image_vector = @(Bitmap) double(reshape(Bitmap,1,row*col));
vector_image = @(Vec) reshape( uint8( min(max(Vec,0),255) ), row, col);
vector_render = @(Vec) imshow(vector_image(Vec));
image_resize = @(Bitmap) [zeros(64, 5) imresize(Bitmap, [64 54]) zeros(64, 5)];
samplesize = 150;
% image showing the difference between your average normal face and the one provided
cd new_faces/faces/
filenames = dir('*.jpg');
file_average_face = 'averagefaceimage_seta.jpeg';
[f0, singular_values, eigen_faces] = eigenface(filenames, 30, 150);
average_a = image_resize(imread(file_average_face));
diff = abs(vector_image(f0) - average_a);
imagesc(diff);
cd ../../
% image showing the difference between your average smiling face and the one provided
samplesize = 150;
cd new_faces/smiling_faces/
filenames = dir('*.jpg');
file_average_face = 'averagefaceimage_setb.jpeg';
[f1, singular_values_smiling, eigen_smiling_faces] = eigenface(filenames, 30, samplesize);
average_b = image_resize(imread(file_average_face));
diff = abs(vector_image(f1) - average_b);
imagesc(diff);
cd ../../
1b: Smiling makes a Difference¶
If your mean normal face is $\overline{\mathbf{f}}_0$, and your mean smiling face is $\overline{\mathbf{f}}_1$, render (using imagesc) the difference $\overline{\mathbf{f}}_0-\overline{\mathbf{f}}_1$ (the average difference between normal faces and smiling faces).
% image showing the difference between the average normal face and average smiling face
smiling_diff = f0 - f1;
imagesc(vector_image(smiling_diff));
1c: Scree Plots and $k$-Approximation¶
Using your downsampled images, perform PCA on each set of faces (normal and smiling).
For each image (normal or smiling), construct its $64^2 \times 1$ vector $\mathbf{f}$. Then, subtract the average face ($\overline{\mathbf{f}}_0$ or $\overline{\mathbf{f}}_1$) and project the remainder on the first $k = 60$ eigenfaces. For example, with a smiling face, the projection of $\mathbf{f}$ on the $j$-th smiling eigenface $\mathbf{e}_j$ is $$ c_j ~=~ {(\mathbf{f} \, - \, \overline{\mathbf{f}}_1)}' \, {\mathbf{e}_j} ~~~~~~~~~~ (j = 1,\dots,k). $$
For each set of faces (normal or smiling), make one large scree plot for the set, showing all sequences of the first $k$ coefficients for each image overplotted (e.g. with hold on).
% (overlaid) scree plots for normal faces
row = 64;
col = 64;
cd new_faces/faces/
filenames = dir('*.jpg');
filesize = size(filenames, 1);
F_norm = zeros(samplesize, row * col); % the array of sample images (stored as vectors)
%numfiles = size(filenames,1);
%rp = randperm(numfiles); % random permutation of the list of image filenames
%sample = rp(1:samplesize); % use the first <samplesize> images as our sample
image_resize = @(Bitmap) [zeros(64, 5) imresize(Bitmap, [64 54]) zeros(64, 5)];
for i = 1 : filesize
Image_File = filenames(i).name;
Face_Matrix = image_resize(imread(Image_File));
F_norm(i,:) = image_vector(Face_Matrix); % the i-th row of F is the i-th image
end
f0 = sum(F_norm,1) / filesize;
%----------------------------------------------------------------------------------
% Subtract the average face from each face in F
%----------------------------------------------------------------------------------
for i = 1 : filesize
F_norm(i,:) = F_norm(i,:) - f0;
end
%----------------------------------------------------------------------------------
% Obtain the eigenfaces U with PCA
%----------------------------------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [U, S, V] = svds( cov(F), k); %%%% expensive
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% singular_values = diag(S);
% do these things more efficiently:
k = 60;
[eigen_faces_60, singular_values_60] = more_efficient_pca( F_norm, k );
%----------------------------------------------------------------------------------
% The columns of the U matrix obtained with PCA are our "eigenfaces".
%
% That is, e_1, e_2, ..., e_k are the first k eigenfaces:
%
% U = ( e_1 | e_2 | ... | e_k | ... )
%
% These eigenfaces are orthonormal, and define our axes of "face space".
%
%----------------------------------------------------------------------------------
%----------------------------------------------------------------------------------
% Finding the eigenface representation of a face
%----------------------------------------------------------------------------------
%
% Any face image (here, a row vector) f can be represented as a sum
%
% f = fbar + c_1 * e_1' + c_2 * e_2' + ... + c_k * e_k'
%
% where the c_j's are the "coordinates in face space" for f, or equivalently,
% if c = [c_1; c_2; ...; c_k; 0; ...; 0] is a row vector of coefficients then
%
% f = fbar + c * U'.
%
% Thus an image f (a 64 x 64 array) is represented by a vector c
% containing only k numbers. Impressive compression!
%
%----------------------------------------------------------------------------------
% t = f - fbar; % remove the average face before finding eigencoefficients
% c = t * U; % eigencoefficients for f: c_j = t * e_j
% c: 1 x k
C_norm = zeros(filesize, k);
figure
for i = 1 : filesize
C_norm(i, :) = F_norm(i,:) * eigen_faces_60;
plot(C_norm(i, 1 : k), '.');
hold on;
end
cd ../../
% (overlaid) scree plots for smiling faces
row = 64;
col = 64;
cd new_faces/smiling_faces/
filenames = dir('*.jpg');
filesize = size(filenames, 1);
F_smiling = zeros(samplesize, row * col); % the array of sample images (stored as vectors)
%numfiles = size(filenames,1);
%rp = randperm(numfiles); % random permutation of the list of image filenames
%sample = rp(1:samplesize); % use the first <samplesize> images as our sample
image_resize = @(Bitmap) [zeros(64, 5) imresize(Bitmap, [64 54]) zeros(64, 5)];
for i = 1 : filesize
Image_File = filenames(i).name;
Face_Matrix = image_resize(imread(Image_File));
F_smiling(i,:) = image_vector(Face_Matrix); % the i-th row of F is the i-th image
end
f1 = sum(F_smiling,1) / filesize;
%----------------------------------------------------------------------------------
% Subtract the average face from each face in F
%----------------------------------------------------------------------------------
for i = 1 : filesize
F_smiling(i,:) = F_smiling(i,:) - f1;
end
%----------------------------------------------------------------------------------
% Obtain the eigenfaces U with PCA
%----------------------------------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [U, S, V] = svds( cov(F), k); %%%% expensive
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% singular_values = diag(S);
% do these things more efficiently:
k = 60;
[eigen_smiling_faces_60, singular_values_smiling_60] = more_efficient_pca( F_smiling, k );
%----------------------------------------------------------------------------------
% The columns of the U matrix obtained with PCA are our "eigenfaces".
%
% That is, e_1, e_2, ..., e_k are the first k eigenfaces:
%
% U = ( e_1 | e_2 | ... | e_k | ... )
%
% These eigenfaces are orthonormal, and define our axes of "face space".
%
%----------------------------------------------------------------------------------
%----------------------------------------------------------------------------------
% Finding the eigenface representation of a face
%----------------------------------------------------------------------------------
%
% Any face image (here, a row vector) f can be represented as a sum
%
% f = fbar + c_1 * e_1' + c_2 * e_2' + ... + c_k * e_k'
%
% where the c_j's are the "coordinates in face space" for f, or equivalently,
% if c = [c_1; c_2; ...; c_k; 0; ...; 0] is a row vector of coefficients then
%
% f = fbar + c * U'.
%
% Thus an image f (a 64 x 64 array) is represented by a vector c
% containing only k numbers. Impressive compression!
%
%----------------------------------------------------------------------------------
% t = f - fbar; % remove the average face before finding eigencoefficients
% c = t * U; % eigencoefficients for f: c_j = t * e_j
% c: 1 x k
C_smiling = zeros(filesize, k);
figure
for i = 1 : filesize
C_smiling(i, :) = F_smiling(i,:) * eigen_smiling_faces_60;
plot(C_smiling(i, :), '.');
hold on;
end
cd ../../
1d: Unusualness of a Face¶
Let us say the unusualness of a face is the $L_2$ norm of its eigenface-coefficient vector (using the first $k=60$ eigenfaces).
Determine, for each set (normal or smiling), the most unusual face.
% the most unusual normal face
%L2Norm = zeros(filesize, 1);
%for i = 1 : samplesize
% L2Norm(i) = norm(C_norm(i, :));
%end
L2Norm = sqrt(sum(C_norm .^ 2, 2));
[maxSum unusualNormalIndex] = max(L2Norm);
% unusualNormalIndex
most_unusual_normal_face = F_norm(unusualNormalIndex, :) + f0;
figure
vector_render(most_unusual_normal_face);
% the most unusual smiling face
%L2Norm = zeros(samplesize, 1);
%for i = 1 : samplesize
% L2Norm(i) = norm(C_smiling(i, :));
%end
L2Norm = sqrt(sum(C_smiling .^ 2, 2));
[maxSum unusualSmilingIndex] = max(L2Norm);
% unusualSmilingIndex
most_unusual_smiling_face = F_smiling(unusualSmilingIndex, :) + f1;
figure
vector_render(most_unusual_smiling_face);
Problem 2: Face Classifiers¶
Develop two different face classifiers using the eigenfaces you've computed; each should be a function that, given a face image $\mathbf{f}$ as input, yields the output value 1 if $\mathbf{f}$ is smiling, and 0 otherwise. (NOTE: or vice-versa; we just need the function to be a classifier)
Specifically, implement the following 3 classifiers that take an input image $\mathbf{f}$:
- {Classifier X}: yield 1 if the normal face unusualness of $\mathbf{f}$ is greater than smiling face unusualness of $\mathbf{f}$, else 0.
- {Classifier Y}: yield 1 if $||{\mathbf{f} - \overline{\mathbf{f}}_0}||^2 ~\geq ~ ||{\mathbf{f} - \overline{\mathbf{f}}_1}||^2$, else 0.
2a: Unusual Face Classification¶
Using each of these classifiers, determine the classification it yields for the two most unusual images you found in the previous question.
% X, Y classifications of the most unusual normal face
classifierX = @(ivec, f0, f1) norm((ivec - f0) * eigen_faces_60) > norm((ivec - f1) * eigen_smiling_faces_60);
classifierX(most_unusual_normal_face, f0, f1)
classifierY = @(ivec, f0, f1) norm((ivec - f0), 1)^2 >= norm((ivec - f1), 1)^2;
classifierY(most_unusual_normal_face, f0, f1)
% X, Y classifications of the most unusual smiling face
classifierX = @(ivec, f0, f1) norm((ivec - f0) * eigen_faces_60) > norm((ivec - f1) * eigen_smiling_faces_60);
classifierX(most_unusual_smiling_face, f0, f1)
classifierY = @(ivec, f0, f1) norm((ivec - f0), 1)^2 >= norm((ivec - f1), 1)^2;
classifierY(most_unusual_smiling_face, f0, f1)
2b: Splitting into Training and Test sets¶
Write a function [Sublist1 Sublist2] = randsplit(List) that takes an array List of length n and splits it randomly into two sublists of size floor(n/2) and ceil(n/2). (Hint: randperm)
Use randsplit to split each of the 200-face sets into a training subset and testing subset of equal size.
For both sets of faces (100 normal faces and 100 smiling faces), compute the average normal and smiling faces $\overline{\mathbf{f}}_0$ and $\overline{\mathbf{f}}_1$ for the training set.
% The randsplit function
function [Sublist1, Sublist2] = randsplit(List)
n = size(List, 1);
IndexArr = randperm(n);
n_floor = int32(floor(n / 2));
Sublist1 = zeros(n_floor, 1);
for i = 1 : n_floor
Sublist1(i) = List(IndexArr(i));
end
n_ceil = int32(ceil(n / 2));
Sublist2 = zeros(n_ceil, 1);
for i = 1 : n_ceil
Sublist2(i) = List(IndexArr(i + n_floor));
end
end
% The average smiling face (for the training set)
samplesize = 200;
k = 60;
cd new_faces/smiling_faces/
filenames = dir('*.jpg');
index = linspace(1, samplesize, samplesize)';
[TrainingIndex, TestingIndex] = randsplit(index);
trainingSize = size(TrainingIndex, 1);
TrainingSmiling = filenames(1 : trainingSize);
for i = 1 : trainingSize
TrainingSmiling(i) = filenames(i);
end
testingSize = size(TestingIndex, 1);
TestingSmiling = filenames(1 : testingSize);
for i = 1 : testingSize
TestingSmiling(i) = filenames(i);
end
[f1bar, singular_values, eigen_smiling_faces_60] = eigenface(TrainingSmiling, k, trainingSize);
figure
vector_render(f1bar);
cd ../../
% The average normal face (for the training set)
samplesize = 200;
k = 60;
cd new_faces/faces/
filenames = dir('*.jpg');
index = linspace(1, samplesize, samplesize)';
[TrainingIndex, TestingIndex] = randsplit(index);
trainingSize = size(TrainingIndex, 1);
TrainingNormal = filenames(1 : trainingSize);
for i = 1 : trainingSize
TrainingNormal(i) = filenames(i);
end
testingSize = size(TestingIndex, 1);
TestingNormal = filenames(1 : testingSize);
for i = 1 : testingSize
TestingNormal(i) = filenames(i);
end
[f0bar, singular_values, eigen_faces_60] = eigenface(TrainingNormal, k, trainingSize);
figure
vector_render(f0bar);
cd ../../
2c: Classifier Error Rate¶
For each of the Classifiers (X, Y), using the average faces you just computed:
- classify each of the 200 faces $\mathbf{f}$ in the testing set, and count classification errors.
- compute the error rate (percentage of errors in test face classifications) for the Classifier.
Then rank the classifiers by their error rate.
For normal faces (using the test set):¶
% X, Y error rates
cd new_faces/faces/
for i = 1 : testingSize
Image_File = TestingNormal(i).name;
Face_Matrix = image_resize(imread(Image_File));
F_norm(i,:) = image_vector(Face_Matrix); % the i-th row of F is the i-th image
end
results = zeros(testingSize, 1);
for i = 1 : testingSize
results(i) = classifierX(F_norm(i, :), f0bar, f1bar);
end
error_rates_X_norm = sum(results) / testingSize
results = zeros(testingSize, 1);
for i = 1 : testingSize
results(i) = classifierY(F_norm(i, :), f0bar, f1bar);
end
error_rates_Y_norm = sum(results) / testingSize
cd ../../
% which of the classifiers has lowest error rate for normal faces in the test set?
% Classifier Y
For smiling faces (using the test set):¶
% X, Y error rates
cd new_faces/smiling_faces/
for i = 1 : testingSize
Image_File = TestingSmiling(i).name;
Face_Matrix = image_resize(imread(Image_File));
F_smiling(i,:) = image_vector(Face_Matrix); % the i-th row of F is the i-th image
end
results = zeros(testingSize, 1);
for i = 1 : testingSize
results(i) = classifierX(F_smiling(i, :), f0bar, f1bar);
end
error_rates_X_norm = sum(~results) / testingSize
results = zeros(testingSize, 1);
for i = 1 : testingSize
results(i) = classifierY(F_smiling(i, :), f0bar, f1bar);
end
error_rates_Y_norm = sum(~results) / testingSize
cd ../../
% which of the classifiers has lowest error rate for smiling faces in the test set?
% Classifier Y
Problem 3: Face Compression¶
In the previous problem you computed the first 60 Eigenface coefficients, and used these to find the most unusual face.
For each $64 \times 64$ image $X$ from your (downsampled) smiling faces, compute the following sequences:
- (descendingly sorted absolute values of) the first 60 Eigenface coefficients for $X$
- (descendingly sorted absolute values of) the first 60 coefficient values from the two-sided FFT of $X$ (in Matlab: fft2(X))
- (descendingly sorted absolute values of) the first 60 coefficient values from the two-sided DCT of $X$ (in Matlab: dct2(X))
- the first 60 singular values from the SVD of $X$.
We get an image compression scheme if we keep only the first $k \leq 60$ coefficients, and discard the rest.
Define $$ \mbox{$k$-coefficient compression error} ~~=~~ \frac{\mbox{(the sum of absolute values of all coefficients after the first $k$)}}{\mbox{(the sum of absolute values of all coefficients)}} . $$
Compute the $k$-coefficient compression error for each of the 4 transforms, $1 \leq k \leq 60$, for the smiling test set.
Rank the 4 transforms above by their average compression error (for $k \leq 60$).
% plot of the Eigenface's k-coefficient compression error (for k <= 60)
%----------------------------------------------------------------------------------
% Obtain the eigenfaces U with PCA
%----------------------------------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [U, S, V] = svds( cov(F), k); %%%% expensive
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% singular_values = diag(S);
% do these things more efficiently:
k = 60;
[eigen_smiling_faces_sample, singular_values_smiling_sample] = more_efficient_pca( F_smiling, testingSize );
C_smiling = zeros(testingSize, testingSize);
for i = 1 : testingSize
C_smiling(i, :) = F_smiling(i,:) * eigen_smiling_faces_sample;
end
k_coeff_comp_error_eigen = zeros(testingSize , k);
for i = 1 : testingSize
for j = 1 : k
% (the sum of absolute values of all coefficients after the first k) /
% (the sum of absolute values of all coefficients)
k_coeff_comp_error_eigen(i, j) = sum(abs(C_smiling(i, j + 1 : testingSize))) / sum(abs(C_smiling(i, 1 : testingSize)));
end
end
figure
plot(1 : k, k_coeff_comp_error_eigen);
title('Eigenfaces compression error')
xlabel('using the first k eigenfaces')
% plot of the two-sided FFT's k-coefficient compression error (for k <= 60)
%%% You might do something like this:
% TwoSidedFFTofX = fft2(X);
% SortedAbsoluteValuesOfFourierCoefficients = sort(abs(TwoSidedFFTofX(:)), 'descend');
% figure
% plot( SortedAbsoluteValuesOfFourierCoefficients(1:60) )
k = 60;
C_smiling = zeros(testingSize, 4096);
for i = 1 : testingSize
C_smiling(i, :) = sort(abs(fft2(F_smiling(i, :))), 'descend');
end
k_coeff_comp_error_fft = zeros(testingSize , k);
for i = 1 : testingSize
for j = 1 : k
% (the sum of absolute values of all coefficients after the first k) /
% (the sum of absolute values of all coefficients)
k_coeff_comp_error_fft(i, j) = sum(abs(C_smiling(i, j + 1 : testingSize))) / sum(abs(C_smiling(i, 1 : testingSize)));
end
end
figure
plot(1 : k, k_coeff_comp_error_fft);
title('FFT compression error')
xlabel('using the first k eigenfaces')
% plot of the two-sided DCT's k-coefficient compression error (for k <= 60)
k = 60;
C_smiling = zeros(testingSize, 4096);
for i = 1 : testingSize
C_smiling(i, :) = sort(abs(dct2(F_smiling(i, :))), 'descend');
end
k_coeff_comp_error_dct = zeros(testingSize , k);
for i = 1 : testingSize
for j = 1 : k
% (the sum of absolute values of all coefficients after the first k) /
% (the sum of absolute values of all coefficients)
k_coeff_comp_error_dct(i, j) = sum(abs(C_smiling(i, j + 1 : testingSize))) / sum(abs(C_smiling(i, 1 : testingSize)));
end
end
figure
plot(1 : k, k_coeff_comp_error_dct);
title('DCT compression error')
xlabel('using the first k eigenfaces')
% plot of the rank-k singular value compression error (for k <= 60)
k = 60;
col = 64;
C_smiling = zeros(testingSize, col);
for i = 1 : testingSize
[U, S, V] = svd(double(vector_image(F_smiling(i,:))));
C_smiling(i, :) = diag(S);
end
k_coeff_comp_error_svd = zeros(testingSize , k);
for i = 1 : testingSize
for j = 1 : k
% (the sum of absolute values of all coefficients after the first k) /
% (the sum of absolute values of all coefficients)
k_coeff_comp_error_svd(i, j) = sum(abs(C_smiling(i, j + 1 : col))) / sum(abs(C_smiling(i, 1 : col)));
end
end
figure
plot(1 : k, k_coeff_comp_error_svd);
title('SVD compression error')
xlabel('using the first k eigenfaces')
% which of the 4 compression schemes has lowest average compression error?
[val index] = min([mean(k_coeff_comp_error_eigen); mean(k_coeff_comp_error_fft); mean(k_coeff_comp_error_dct); mean(k_coeff_comp_error_svd)]);
if index == 1
fprintf('Eigenfaces');
elseif index == 2
fprintf('FFT');
elseif index == 3
fprintf('DCT');
else
fprintf('SVD');
end
fprintf(' has the smallest average compression error');