Load in needed libraries

library(jpeg)
library(foreach)
library(EBImage)
# Set files to be all the shoe images
files=list.files(path='jpg/',pattern="\\.jpg", full.names=TRUE)
# print the dimenstions to get the height and width
print(dim(readJPEG(files[1])))
## [1] 1200 2500    3
# Set the height, width and scale variables
height=1200
width=2500
scale=20
# Loop over the images plot, and resize them and add them to the rezized images list

resized_images <- list()

# Initialize a list to store flattened image vectors
image_vectors <- list()
par(mfrow=c(6,3)) #set graphics to 6 x 3 table
par(mai=c(.3,.3,.3,.3)) #set margins 
# Loop through image files, flatten each image, and store in the list
for (image_file in files) {
  # Read the image
  img <- readJPEG(image_file)
  res = dim(img)[2:1]
  # Get the dimensions of the image
  dimensions <- dim(img)
  width <- dimensions[2]
  height <- dimensions[1]
  
  # Resize the image
  # Flatten the image into a vector and store in the list
  image_vectors[[image_file]] <- array(EBImage::resize(img,height/scale, width/scale),dim=c(1, height/scale, width/scale,3))
  
  plot(1,1,xlim=c(1,res[1]),ylim=c(1,res[2]), #set the X Limits by size
         asp=1, #aspect ratio
         type='n', #don't plot
         xaxs='i',yaxs='i',#prevents expanding axis windows +6% as normal
         xaxt='n',yaxt='n',xlab='',ylab='', # no axes or labels
         bty='n') # no box around graph
  rasterImage(img,1,1,res[1],res[2]) #image, xleft,ybottom,xright,ytop
}

# Stack resized image vectors into a matrix
image_matrix <- do.call(rbind, image_vectors)
# Setup the newdata matrix for pca, run and get principal components
newdata=image_matrix
dim(newdata)=c(length(files),height*width*3/scale^2)
mypca=princomp(t(as.matrix(newdata)), scores=TRUE, cor=TRUE)
# Get the princpal componets and save them to a variable
mycomponents=mypca$sdev^2/sum(mypca$sdev^2)
mycomponents
##      Comp.1      Comp.2      Comp.3      Comp.4      Comp.5      Comp.6 
## 0.692820228 0.101224662 0.051062366 0.027277491 0.018999399 0.016249609 
##      Comp.7      Comp.8      Comp.9     Comp.10     Comp.11     Comp.12 
## 0.013994442 0.012060722 0.009698242 0.009058322 0.008458196 0.007987044 
##     Comp.13     Comp.14     Comp.15     Comp.16     Comp.17 
## 0.007632793 0.006697401 0.006182362 0.005955453 0.004641268
# Set threshold to 80% variance
# loop over the compoents and sum them until we reach 80% variance then print how many components were used
threshold <- .80
cumulative_variance <- 0
num_components <- 0
for(comp in mycomponents) {
  cumulative_variance <- cumulative_variance + comp
  num_components <- num_components +1
  if(cumulative_variance >= threshold) {
    print(num_components)
    break;
  }
}
## [1] 3
#. Show the variance for first 2 components
sum(mycomponents[1:2]) 
## [1] 0.7940449
#. Show variance for first 3 components
sum(mycomponents[1:3]) 
## [1] 0.8451073
# Now plot the 17 eigen shoes using the principal components
mypca2=t(mypca$scores)
dim(mypca2)=c(length(files),height/scale,width/scale,3)
par(mfrow=c(5,5))
par(mai=c(.001,.001,.001,.001))
for (i in 1:17){  #plot the first 81 Eigenshoes only
  path = writeJPEG(mypca2[i,,,], quality=1,bg="white")
  jpg = readJPEG(path, native=T) # read the file
  res = dim(jpg)[2:1] # get the resolution, [x is 2, y is 1]
  # initialize an empty plot area if add==FALSE
  plot(1,1,xlim=c(1,res[1]),ylim=c(1,res[2]), #set the X Limits by size
       asp=1, #aspect ratio
       type='n', #don't plot
       xaxs='i',yaxs='i',#prevents expanding axis windows +6% as normal
       xaxt='n',yaxt='n',xlab='',ylab='', # no axes or labels
       bty='n') # no box around graph
  rasterImage(jpg,1,1,res[1],res[2]) #image, xleft,ybottom,xright,ytop
}