library(doParallel)
## Loading required package: foreach
## Loading required package: iterators
## Loading required package: parallel
library(foreach)
library(jpeg)
library(EBImage)
library(OpenImageR)
## 
## Attaching package: 'OpenImageR'
## The following objects are masked from 'package:EBImage':
## 
##     readImage, writeImage
library(utils)

Load images and set the height, width and scale

shoes = list.files(pattern="\\.jpg$",full.names = TRUE)
num_files = length(shoes)
height <- 1200
width <- 2500
scale <- 20

Create an empty array to load the images and resize to the scaled dimensions

im=array(rep(0,length(shoes)*height/scale*width/scale*3), dim=c(length(shoes), height/scale, width/scale,3))

for (i in 1:length(shoes)){
  temp=resize(readJPEG(shoes[i]),height/scale, width/scale)
  im[i,,,]=array(temp,dim=c(1, height/scale, width/scale,3))}

Function to view and plot the images by looping through each image in the matrix

plot_jpeg = function(path, add=FALSE)
{ jpg = readJPEG(path, native=T) # read the file
  res = dim(jpg)[2:1] # get the resolution, [x, y]
  if (!add) # initialize an empty plot area if add==FALSE
    plot(1,1,xlim=c(1,res[1]),ylim=c(1,res[2]),asp=1,type='n',xaxs='i',yaxs='i',xaxt='n',yaxt='n',xlab='',ylab='',bty='n')
  rasterImage(jpg,1,1,res[1],res[2])
}

par(mfrow=c(3,6))
par(mai=c(.01,.01,.01,.01))
for (i in 1:length(shoes)){  #plot the first images only
  plot_jpeg(writeJPEG(im[i,,,]))
}

## Load the transposed image RGB vectors into a flat data frame for linear transformation

flat <- matrix(0, length(shoes), prod(dim(im))) 
for (i in 1:length(shoes)) {
  newim <- readJPEG(shoes[i])
  r <- as.vector(im[i,,,1])
  g <- as.vector(im[i,,,2])
  b <- as.vector(im[i,,,3])
  flat[i,] <- t(c(r, g, b))
}
shoes_df <- as.data.frame(t(flat))

Scale images, create a covariance matrix to determine correlation between each image and determine eigenvalues

scaled_img = scale(shoes_df, center = TRUE, scale = TRUE)
correl = cor(scaled_img)
eigenv = eigen(correl)

Calculate the variance to determine which image accounts for 80% of the variability

var <- cumsum(eigenv$values) / sum(eigenv$values)
var
##  [1] 0.6928202 0.7940449 0.8451073 0.8723847 0.8913841 0.9076338 0.9216282
##  [8] 0.9336889 0.9433872 0.9524455 0.9609037 0.9688907 0.9765235 0.9832209
## [15] 0.9894033 0.9953587 1.0000000
scaling=diag(eigenv$values[1:2]^(-1/2)) / (sqrt(nrow(scaled_img)-1))
eigenshoes=scaled_img%*%eigenv$vectors[,1:2]%*%scaling

The second image accounts for 80% of the variability.

imageShow(array(eigenshoes[,2], c(60,125,3)))