url = "C:/Users/okereann/Desktop/shoes/"
# list the files
img_files <- list.files(url, pattern = "\\jpg$")
# specify dimensions
img_lengths = length(img_files)
img_height=120; 
img_width=250;
fhi = matrix(1, nrow = 3, ncol = 3)
fhi[2, 2] = -8
# Create the plot_jpeg function
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])
}

im=array(rep(0,length(img_files)*img_height*img_width*3), dim=c(img_lengths, img_height, img_width, 3))

# resize the images
for (i in 1:length(img_files)){

  temp= EBImage::resize(readJPEG(paste0(url, img_files[i])), img_height, img_width)
  im[i,,,]=array(temp,dim=c(1, img_height, img_width, 3))
}
flat=matrix(0, img_lengths, prod(dim(im)))
for (i in 1:img_lengths){
  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))
# plot the images
par(mfrow = c(3,3))
par(mai = c(.3,.3,.3,.3))
for (i in 1:img_lengths){ 
plot_jpeg(writeJPEG(im[i,,,]))
}

# Calculate covariance
scaled_shoes_df = scale(shoes_df, center = TRUE, scale = TRUE)
Sigma_= cor(scaled_shoes_df)

# Get the eigen components
eigen_imagery = eigen(Sigma_)
cumsum(eigen_imagery$values) / sum(eigen(Sigma_)$values)
##  [1] 0.6916309 0.7911549 0.8442585 0.8714110 0.8903922 0.9070168 0.9213980
##  [8] 0.9335995 0.9433646 0.9523060 0.9606716 0.9685299 0.9760625 0.9826631
## [15] 0.9891121 0.9950668 1.0000000
scaling = diag(eigen_imagery$values[1:5]^(-1/2)) / (sqrt(nrow(scaled_shoes_df)-1))
eigenshoes = scaled_shoes_df %*% eigen_imagery$vectors[,1:5] %*% scaling


# Compute the principal components
new_data = im
dim(new_data) = c(img_lengths, img_height * img_width*3)
pca_of_new_data = princomp(t(as.matrix(new_data)), scores=TRUE, cor=TRUE)

# get pca scores
pca_scores = t(pca_of_new_data$scores)
dim(pca_scores) = c(img_lengths, img_height, img_width, 3)


# plot the new images
par(mfrow=c(5,5))
par(mai=c(.001,.001,.001,.001))
for (i in 1:img_lengths){
plot_jpeg(writeJPEG(pca_scores[i,,,], bg="white"))  
}