Lesson 10.2 - PCA for Facial Recognition

Load Packages

```r
library(caret)
library(shape)

<!-- rnb-source-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->




<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxubGlicmFyeShoZXhWaWV3KVxucmF3X2RhdGEgPC0gcmVhZFJhdyhcImRhdGEvbWl0X2ZhY2VzL3Jhd2RhdGEvMjIyN1wiKVxuaW1hZ2VfbWF0cml4IDwtIG1hdHJpeChyYXdfZGF0YSRmaWxlUmF3LCBucm93ID0gMTI4LCBcbiAgICAgICAgICAgICAgICAgICAgICAgbmNvbCA9IDEyOCwgYnlyb3cgPSBUUlVFKVxucGxvdChhcy5yYXN0ZXIoaW1hZ2VfbWF0cml4KSlcbmBgYCJ9 -->

```r
library(hexView)
raw_data <- readRaw("data/mit_faces/rawdata/2227")
image_matrix <- matrix(raw_data$fileRaw, nrow = 128, 
                       ncol = 128, byrow = TRUE)
plot(as.raster(image_matrix))

set.seed(1)
par(mfrow=c(10,10))
for (i in sample(1:3991, 100)){
  par(mar=c(0,0,0,0))
  image_matrix <- matrix(faces[i,], nrow = 128, 
                         ncol = 128, byrow = TRUE)
  plot(as.raster(image_matrix))
}
par(mfrow=c(1,1))

temp_df <- data.frame(faces)
df_faces <- temp_df
for (i in 1:ncol(temp_df)){
  df_faces[,i] <- as.numeric(temp_df[,i])
}

face_pca <- prcomp(df_faces)
dim(face_pca$x)

[1] 3991 3991

dim(face_pca$rotation)

[1] 16384  3991

ev <- face_pca$sdev**2 / sum(face_pca$sdev**2)
cum_ev <- cumsum(ev)

plot(1:1000, cum_ev[1:1000], cex=0.25, pch='',
     xlab='Number of Principal Components', 
     ylab='Cumulative Variance Explained')
lines(1:1000, cum_ev[1:1000], col='steelblue', lwd=2)

abline(v=200, lty=3, lwd=2, col='steelblue')
abline(h=cum_ev[200], lty=3, lwd=2, col='steelblue')

set.seed(1)
par(mfrow=c(5,5))
for (i in 1:25){
  par(mar=c(0,0,0,0))
  image_matrix <- matrix(face_pca$rotation[,i], nrow = 128, 
                           ncol = 128, byrow = TRUE)
  
  minval <- min(image_matrix)
  maxval <- max(image_matrix)
  
  image_matrix <- (image_matrix - minval) / (maxval - minval)
  plot(as.raster(image_matrix))
}
par(mfrow=c(1,1))

set.seed(1)
par(mfrow=c(5,5))
for (i in 26:50){
  par(mar=c(0,0,0,0))
  image_matrix <- matrix(face_pca$rotation[,i], nrow = 128, 
                           ncol = 128, byrow = TRUE)
  
  minval <- min(image_matrix)
  maxval <- max(image_matrix)
  
  image_matrix <- (image_matrix - minval) / (maxval - minval)
  plot(as.raster(image_matrix))
}
par(mfrow=c(1,1))

ix <- sample(1:3991, 1)
ix <- 2177
par(mar=c(0,0,0,0))
image_matrix <- matrix(faces[ix,], nrow = 128, 
                       ncol = 128, byrow = TRUE)
plot(as.raster(image_matrix))

face_pca$x[2177,1:25]

        PC1         PC2         PC3         PC4         PC5         PC6         PC7 
-3345.07194   111.59096  1145.27932 -1413.52145  1812.93469 -2120.45176   528.41917 
        PC8         PC9        PC10        PC11        PC12        PC13        PC14 
  356.65318  1097.28055  -847.17761    88.20587    36.28544   694.17292   -75.62807 
       PC15        PC16        PC17        PC18        PC19        PC20        PC21 
 -576.03196  -544.69454  -117.88797   320.06718  -227.21212  -526.35007   353.70063 
       PC22        PC23        PC24        PC25 
  492.15796  -278.60709  -239.30506    72.74486 

cur_img <- face_pca$rotation[,1] * 0

par(mfrow=c(40,5))
for (i in 1:200){
  par(mar=c(0,0,0,0))
  
  cur_img <- cur_img + face_pca$x[ix, i] * face_pca$rotation[,i]
  
  image_matrix <- matrix(cur_img, nrow = 128, 
                           ncol = 128, byrow = TRUE)
  
  minval <- min(image_matrix)
  maxval <- max(image_matrix)
  image_matrix <- (image_matrix - minval) / (maxval - minval)
  plot(as.raster(image_matrix))
}
par(mfrow=c(1,1))

ix <- sample(1:3991, 1)
#ix <- 2177
par(mar=c(0,0,0,0))
image_matrix <- matrix(faces[ix,], nrow = 128, 
                       ncol = 128, byrow = TRUE)
plot(as.raster(image_matrix))

face_pca$x[2177,1:25]

cur_img <- face_pca$rotation[,1] * 0

par(mfrow=c(40,5))
for (i in 1:200){
  par(mar=c(0,0,0,0))
  
  cur_img <- cur_img + face_pca$x[ix, i] * face_pca$rotation[,i]
  
  image_matrix <- matrix(cur_img, nrow = 128, 
                           ncol = 128, byrow = TRUE)
  
  minval <- min(image_matrix)
  maxval <- max(image_matrix)
  image_matrix <- (image_matrix - minval) / (maxval - minval)
  plot(as.raster(image_matrix))
}
par(mfrow=c(1,1))