Session 7: Eigenfaces

# These are the standard imports for CS 111
# This list may change as the quarter goes on

import sys
sys.path.append("/Users/youngmin/Documents/CS/Introduction to Computational Science - F21/CS111-2021-fall/Python")

import os
import time
import math
import numpy as np
import numpy.linalg as npla
from scipy import linalg as spla
import scipy.sparse
import scipy.sparse.linalg
from scipy import integrate
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import axes3d
%matplotlib inline
import cs111

Faces as Matrices

# Loading a black-and-white face image for its visualization.
Face = plt.imread( "faces/face000.bmp", format="bmp" )
print( "Face image has shape:", Face.shape )
print( "            and type:", Face.dtype )      # Face is a 64x64 matrix of gray scale
print( "       Minimum value:", np.min( Face ) )  # values thay may range between 0 and 255.
print( "       Maximum value:", np.max( Face ) )

# Let's show the image.
plt.figure( figsize=(3,3) )
plt.title( r"A $64\times 64$ black-and-white face image" )
plt.imshow( Face, cmap="gray" )

Face image has shape: (64, 64)
            and type: uint8
       Minimum value: 0
       Maximum value: 182

<matplotlib.image.AxesImage at 0x123c78d90>

Representing a Data Set of Faces as a Matrix

We will create a data set $F$ of gray-scale images. Each image will be a row in $F$, and their pixels will be placed in the columns.

The accompanying folder faces/ has a data set of 178 black-and-white face images. Out of these 178 faces, there are a few of them (indicated in face_descriptions.py) that should be omitted. Each image is of size $64\times 64$ pixels.

To build our face data set matrix $F$,

1. Reshape each $i^{th}$ image into a row vector, $\mathbf{f}_i^T$, and
2. Stack all of our $1\times 64^2$ vectors in $F$.

First, let's define some utility functions. These functions assume that images are gray-scale matrices of uint8 values between 0 and 255, and vectors are double precision so that we can perform arithmetics on them.

# Some global variables.
HEIGHT = WIDTH = 64

################# Useful functions to go from images to vectors and back #################
imageToVector = lambda Bitmap: Bitmap.reshape( -1 ).astype( np.float64 )
vectorToImage = lambda Vector: Vector.clip( 0, 255 ).reshape( HEIGHT, WIDTH ).astype( np.uint8 )
renderVector = lambda Vector: plt.imshow( vectorToImage( Vector ), cmap="gray" )
scaleImageIntensities = lambda A: np.round( (A - np.min( A )) / (np.max( A ) - np.min( A )) * 255 ).astype( np.uint8 )
##########################################################################################

Building Matrix $F$

By using these functions we can build the $n$-by-$p$ matrix $F$ with $n = 169$ faces stacked as row vectors with $p = 64^2$ columns. This way, the image vector $\mathbf{f}_i^T$ is at the $i^{th}$ row of $F$:

$$ F = \begin{array}{r} \textrm{face image 0} \\ \textrm{face image 1} \\ \vdots \\ \textrm{face image }n - 1 \\ \end{array}\begin{pmatrix} f_{0,0} & f_{0,1} & \cdots & f_{0,p-1} \\ f_{1,0} & f_{1,1} & \cdots & f_{1,p-1} \\ \vdots & \vdots & \ddots & \vdots \\ f_{n-1,0} & f_{n-1,1} & \cdots & f_{n-1,p-1} \\ \end{pmatrix} \in \mathbb{R}^{169\times 64^2}. $$

Let's next load all images in the faces/ directory to build $F$.

import face_descriptions   # Includes the file names we are interested in.

F = []                     # We append to F first as a list; then we make it into a matrix.
for fileName in sorted( face_descriptions.face_features ):
    if fileName not in face_descriptions.image_to_omit:
        bitMap = plt.imread( "faces/" + fileName )
        F.append( imageToVector( bitMap ) )        # Flatenning the image.
F = np.array( F )          # Now F is a matrix.

print( "Shape of F is:", F.shape )

plt.figure( figsize=(10, 3) )
plt.title( r"Entire matrix $F$ with $169\times 64^2$ pixels in it" )
plt.imshow( F, cmap="gray" )

Shape of F is: (169, 4096)

<matplotlib.image.AxesImage at 0x123d5bf70>

The Average Face

We can compute the average face as $\mathbf{m} = \frac{1}{n}\sum_{i}\mathbf{f}_i$.

m = np.mean( F, axis=0 )   # Notice the axis=0 parameter!

# Plotting resulting mean face.
plt.figure( figsize=(3,3) )
plt.title( r"The mean face" )
renderVector( m )

<matplotlib.image.AxesImage at 0x123ec8eb0>

PCA of the Face Matrix

Let $\mathbf{m}^T$ be the mean row vector of matrix $F$ as defined above. Suppose we stack $\mathbf{m}^T$, yet again, $n$ times in a matrix $M$. Then, the covariance matrix of $F$ is given by

$$C =\textrm{cov}(F) = \frac{1}{n-1}(F - M)^T (F - M)$$

If we define $X = (F - M)$, then

$$\textrm{cov}(X) = \frac{1}{n-1}X^T X = \frac{1}{n-1}(F - M)^T (F - M) = \textrm{cov}(F) = C,$$

and

$$C = USV^T = VSV^T,$$

since $C$ is symmetric positive semidefinite and $U = V$. We call the columns in $V$ the principal components of $F$, and we can use them to express any face $\mathbf{f}$. Suppose $\mathbf{x} = \mathbf{f} - \mathbf{m}$ is a vectorized/flattened face image with all pixel values centered around the corresponding means. We can express $\mathbf{x}$ as

$$\mathbf{x} = V\mathbf{c} = c_0\mathbf{e}_0 + c_1\mathbf{e}_1 + \cdots + c_{p-1}\mathbf{e}_{p-1},$$

where $\mathbf{c} = (c_0, c_1, ..., c_{p-1})^T$ is a vector of coefficients, and

$$ V = \left(\begin{array}{c|c|c|c} ~ & ~ & ~ & ~ \\ \mathbf{e}_0 & \mathbf{e}_1 & \cdots & \mathbf{e}_{p-1} \\ ~ & ~ & ~ & ~ \\ \end{array}\right) $$

is the $p$-by-$p$ matrix of singular vectors of $\textrm{cov}(X)$.

To find the vector of coefficients/projections $\mathbf{c}$, we use

$$\mathbf{c} = V^T \mathbf{x},$$

and so $\mathbf{x} = V\mathbf{c} = V\left(V^T \mathbf{x}\right) = \mathbf{x}$. Furthermore, since $\mathbf{x} = \mathbf{f} - \mathbf{m}$, then

$$\mathbf{x} = (\mathbf{f} - \mathbf{m}) = V\mathbf{c},$$

which implies

$$\mathbf{f} = \mathbf{m} + V\mathbf{c}.$$

This is our cookbook to fabricate faces.

Eigenfaces

With PCA, we can perform data compression or low-rank approximation if we use only the first $k$ singular vectors in $V$. If we do so with any face $\mathbf{f}$, then

$$\mathbf{f}^{(k)} = \mathbf{m} + V^{(k)}\mathbf{c} = \mathbf{m} + c_0\mathbf{e}_0 + c_1\mathbf{e}_1 + \cdots + c_{k-1}\mathbf{e}_{k-1}.$$

The singular vectors in $V$ are called eigenfaces.

# We already know F and m.  Let's compute X: the centered matrix F around the column means.
X = F - m    # Notice that m will be broadcasted to fit the dimensions of X.

# Compute the covariance matrix of the the centered data.
C = np.cov( X.T )
print( "Shape of cov(X):", C.shape )

# Perform PCA on C using the SVD.
U, sigma, Vt = spla.svd( C )
Eigenfaces = Vt.T         # The Eigenfaces are the eigenvectors (U = V).

# Let's look at the first largest k eigen/singular values.
k = 30

# Plotting eigenfaces and the explained variance for first k singular values.
fig = plt.figure( figsize=(10,7) )
plt.subplot( 231 )
plt.title( "Singular values" )
plt.plot( sigma[:k], "." )
plt.subplot( 232 )
plt.title( "First Eigenface" )
renderVector( scaleImageIntensities( Eigenfaces[:,0] ) )
plt.subplot( 233 )
plt.title( "Second Eigenface" )
renderVector( scaleImageIntensities( Eigenfaces[:,1] ) )
plt.subplot( 234 )
plt.title( "Third Eigenface" )
renderVector( scaleImageIntensities( Eigenfaces[:,2] ) )
plt.subplot( 235 )
plt.title( "Fourth Eigenface" )
renderVector( scaleImageIntensities( Eigenfaces[:,3] ) )
plt.subplot( 236 )
plt.title( "Fifth Eigenface" )
renderVector( scaleImageIntensities( Eigenfaces[:,4] ) )
plt.show()

Shape of cov(X): (4096, 4096)

By modifying the coefficient of an eigenface (i.e., $c_j$ for $\mathbf{e}_j$ above), we can alter the appearance of the face. Let's modify the coefficient for the fourth (i.e., $j = 3$) eigenface, which emphasizes the tip of the nose and the eyebrows. In general, the brighter the pixels in an eigenface, the more it emphasizes those pixels/features in a face.

# At this point, we know the mean face m and the Eigenfaces V.
f = plt.imread( "faces/face000.bmp" )    # Let's modify the first face in the data set.
f = imageToVector( f )
x = f - m
c = Vt @ x                               # Equivalently, Eigenfaces.T @ x: How many coefficients are in c?

factors = [1, 3, 6, 9]
plt.figure( figsize=(10,4) )
subplotCode = 140
for i in range( len( factors ) ):
    factor = factors[i]
    plt.subplot( subplotCode + i + 1 )
    cNew = c.copy()
    cNew[3] *= factor
    renderVector( m + Eigenfaces @ cNew ) # Recall that Eigenfaces = V.
    plt.xlabel( "factor = {}".format( factor ) )
    if i == 1:
        plt.title( "A face image with 4th Eigenface coefficient multiplied by a constant factor" )

The $k^{th}$ Approximation of a Face Using Eigenfaces

By having our face cookbook ($V$ = eigenfaces) at hand, we can approximate a face $\mathbf{f}$ with just $k$ numbers!

$$\mathbf{f}^{(k)} = \mathbf{m} + V^{(k)}\mathbf{c} = \mathbf{m} + c_0\mathbf{e}_0 + c_1\mathbf{e}_1 + \cdots + c_{k-1}\mathbf{e}_{k-1}.$$

############## The kth approximation of a face from the cookbook ##############
k = 50                                   # You must use k <= 169 because min(row,col)=169.
f = plt.imread( "faces/face000.bmp" )    # Let's take again the first face in the data set.
f = imageToVector( f )
x = f - m
c = Eigenfaces[:,:k].T @ x

# Plotting the resulting kth approximation and the original image.
plt.figure( figsize=(6,3))
plt.subplot( 121 )
renderVector( f )
plt.xlabel( "Original face" )
plt.subplot( 122 )
renderVector( m + Eigenfaces[:,:k] @ c )
plt.xlabel( "Approximated face" )
plt.suptitle( "Approximating a face using k = {} Eigenfaces".format( k ) )

Text(0.5, 0.98, 'Approximating a face using k = 50 Eigenfaces')

############## The kth approximation of a face *not* using Eigenfaces ##############
k = 20
f = plt.imread( "faces/face000.bmp" )    # Let's take again the first face in the data set.
U2, sigma2, V2t = spla.svd( f )
fk = np.zeros( f.shape )
for i in range( k ):
    fk += sigma2[i] * np.outer( U2[:,i], V2t[i,:] )

plt.figure( figsize=(6,3))
plt.subplot( 121 )
plt.imshow( f, cmap="gray" )
plt.xlabel( "Original face" )
plt.subplot( 122 )
plt.imshow( fk, cmap="gray" )
plt.xlabel( "Approximated face" )
plt.suptitle( "Approximating a face using k = {} singular vector and no eigenfaces".format( k ) )

Text(0.5, 0.98, 'Approximating a face using k = 20 singular vector and no eigenfaces')

Generating Random Faces

By knowing the mean face $\mathbf{m}$ and the eigenfaces in $V$, we can synthesize any face $\tilde{\mathbf{f}}$ by sampling a vector of coefficients $\tilde{\mathbf{c}}$ from some normal distribution:

$$\tilde{\mathbf{f}}^{(k)} = \mathbf{m} + V^{(k)}\tilde{\mathbf{c}} = \mathbf{m} + \tilde{c}_0\mathbf{e}_0 + \tilde{c}_1\mathbf{e}_1 + \cdots + \tilde{c}_{k-1}\mathbf{e}_{k-1}.$$

############## Generating a random face from the cookbook ##############
k = 20                                   # Use only the top k Eigenfaces.
Coeffs = X @ Eigenfaces[:,:k]            # An nxk matrix of face coefficients.
mu = np.mean( Coeffs, axis=0 )           # Mean row vector of coefficients.
std = np.std( Coeffs, ddof=1, axis=0 )   # Standard deviations along each column of coefficients.

print( "Variance from function:" )       # Checking that the variances on coefficients should be the same
print( std ** 2 )                        # as the singular values!
print( "\nVariance from SVD:" )
print( sigma[:k] )

# To generate the new face, we must generate new coefficients based on the
# statistics from the population: the mean and std obtained above.
newFaceCoeffs = std * np.random.randn( k ) + mu  # This creates coefficients with a normal distribution
newFace = m + Eigenfaces[:,:k] @ newFaceCoeffs   # whose mean is mu and variance is sigma^2.

# Plotting new random face.
plt.figure( figsize=(3,3) )
plt.title( "Random face" )
renderVector( newFace )

Variance from function:
[594316.78282257 388739.27150947 220110.71998256 215765.40290314
 148366.11611318 136747.66077291  83176.09921457  70937.49110233
  59014.92801102  55035.63320068  49892.09229926  45004.20836223
  40254.16226681  35670.46406856  34394.65172355  32401.94970633
  30830.37818509  28070.31809972  26477.49452898  24139.03654889]

Variance from SVD:
[594316.78282257 388739.27150947 220110.71998256 215765.40290314
 148366.11611318 136747.66077291  83176.09921457  70937.49110233
  59014.92801102  55035.63320068  49892.09229926  45004.20836223
  40254.16226681  35670.46406856  34394.65172355  32401.94970633
  30830.37818509  28070.31809972  26477.49452898  24139.03654889]

<matplotlib.image.AxesImage at 0x11064daf0>