Introduction

This assignment asks us to build and visualize eigenimagery that accounts for 80% of the variability in a set of JPEG files showing various shoes, all positioned similarly. The task…

Vectorization

Linear Algebra does not deal in images, but in vectors. How, then, can it be used to solve such a visual problem?

Fortunately, digital images can be represented as vectors, where each element represents a degree of colorization or gradation for a particular pixel. Even more fortunately, R contains tools to do exactly that! Once we have our images in vector form, we can act on them with matrix operations just as we would with other vector-based problems. So, let’s begin by building a matrix representing our JPEG images.

First, we can import the jpeg R library, which allows us to import images. Every other function we need can be called natively in R.

library(jpeg)

To test that I can a.) import JPEG images into R and b.) access the folder where my images reside, I will run a test to pull in the first file from list. Note that in the following line of code I am:

  • Using list.files() to pull a list of file names in the jpg folder, then indexing it using [1] to pull only one file name.
  • Pasting ‘jpg’ and the file name together using a ‘/’ separator to properly access the folder and file.
  • Creating a JPEG object in R with the readJPEG function.
my_jpg <- readJPEG(paste('jpg',list.files('jpg')[1],sep='/'))

dim(my_jpg)
## [1] 1200 2500    3

Note that my JPEG has 3 dimensions, reflecting the height, width, and color profile (RGB) of my image respectively.

For the sake of simplicity, I will b working with my images in grayscale. That means I need a repeatable function that can make operate on my image vectors accordingly. The below function multiplies each red, green and blue value for every pixel in the image by certain scalars, resulting in a gray but otherwise recognizable version of an image. This is called the luminosity method and the scalars are constant regardless of the image.

For efficiency, I will also use the same function to “flatten” the image matrix into a single vector. This is important for later steps, which require a single matrix where each row is one vector representing a single image.

turn_img_to_grayscale_vec <- function(img) {
  #these scalars, when multiplied by the R, G and B values, produce a grayscale copy of the image matrix
  grayscale_vec <- 0.299 * img[,,1] + 0.587 * img[,,2] + 0.114 * img[,,3]
  #we return the object as a vector using c()
  return(c(grayscale_vec))
}

now for every image…make the matrix

#create list of file names from jpg directory
file_names <- c(dir('jpg'))

#initialize the matrix with the first filename in the list
img_mtx <- turn_img_to_grayscale_vec(readJPEG(paste('jpg',file_names[1],sep='/')))

#loop through remaining file names, producing flat grayscale vectors for each and binding them
#iteratively to my overall matrix.
for (file_name in tail(file_names, -1)) {
  new_img_vec <- turn_img_to_grayscale_vec(readJPEG(paste('jpg',file_name,sep='/')))
  img_mtx <- rbind(img_mtx,new_img_vec)
}

#display matrix dimensions
dim(img_mtx)
## [1]      17 3000000

We now have a numerical matrix 17 rows long and 3 million columns wide. That makes intuitive sense, since we have 17 images and each image is labeled as 2500 x 1200 pixels, or 3 million total pixels. Each row now represents one of those images (in grayscale), and each column represents a precise pixel location for those images.

The next step is running the matrix through a function that conducts Principal Component Analysis, or PCA.

pca_result <- prcomp(img_mtx, center = TRUE, scale. = FALSE)

What this function has done is essentially “normalize” each row based on the mean pixel values for all images, then reduce the matrix greatly in size to a set of “principal components.” Theoretically, these principal components can account for a high percentage of the variability among our set of images.

We can now find the cumulative variance by dividing the cumulative sum of the square of standard deviations (an iterative way to collect the variances of each principal component) by the overall sum of those squares (the total variance). This gives us the percentage of variance that is explained by each component on step-wise bases.

cumulative_variance <- cumsum(pca_result$sdev^2) / sum(pca_result$sdev^2)

Since we now have an ordered vector representing the cumulative explained variance, we can limit ourselves to the values that get us to 80% explainability.

cumulative_variance
##  [1] 0.4742360 0.5969741 0.6626158 0.7180658 0.7649976 0.8002861 0.8326499
##  [8] 0.8617583 0.8866735 0.9100493 0.9301857 0.9483910 0.9657372 0.9809233
## [15] 0.9939490 1.0000000 1.0000000

We can see from our vector that we reach 80% at our 6th value. Therefore, the “eigenimages” are the first 6 vectors of the rotation matrix from our PCA results:

eigenimages <- pca_result$rotation[, 1:6]