prcomp_mystery.R

library(pixmap)
library(ppls)

# use the working directory of where-ever you are: WD. Set this yourself
# Ex. for me: WD <- "/Users/artsyinc/Downloads/w4240"
data <- read.csv(paste(WD, "/image_matrix2.csv", sep = ""), header = F)
image.matrix <- matrix(unlist(data), nrow = nrow(data), ncol = ncol(data))
# the matrix is features x trials
# each feature is a pixel in an image, each feature is an image of a person
mean.face.v <- as.vector(apply(image.matrix, 2, mean))
image.m.centered <- t(apply(image.matrix, 1, function(row) { row - mean.face.v }))
X <- t(image.m.centered)

########### some necessary methods to setup for later

ith.eigenface.for.image.v <- function(i, v, e.faces) {
  v.centered <- as.vector(v - mean.face.v)
  coef <- v.centered %*% as.vector(e.faces[, i])
  coef * as.vector(e.faces[, i])
}
# generate cumsum of projections onto eigenfaces
cumsum.eigenfaces <- function(top.eigenfaces) {
  top.eigenfaces.cum <- top.eigenfaces # copy
  for (i in seq_along(top.eigenfaces.cum)) {
    sum <- mean.face.v
    for (j in c(1:i)) {
      sum <- sum + top.eigenfaces[[j]]
    }
    top.eigenfaces.cum[[i]] <- sum
  }
  top.eigenfaces.cum
}
# plotting a row image
plot.row.v <- function(row.v, str = "") {
  m <- matrix(row.v, nrow = 192) # all the images are this size
  pix <- pixmapGrey(m)
  plot(pix, main = str)
}

########### prcomp way:
pca <- prcomp(X, center = F)
pca.x.norm <- apply(pca$x, 2, normalize.vector)

########### do it yourself way:
XTX <- t(X) %*% X # missing the 1 / n - 1 for cov matrix b/c we normalize later anyway
eigen <- eigen(XTX)
eigenvectors.XTX.col <- eigen$vectors
# this is equivalent to:
eigenvectors.XTX.col.copy <- -1 * svd(X)$v

principal.component.scores <- apply(eigenvectors.XTX.col, 2, function(c) {
  normalize.vector(X %*% matrix(c, ncol = 1))
}) # same as: X %*% t(eigenvectors.XTX.col) but with normalization

##### test either method's principle component scores:
# run one of these lines:
# eigenfaces.as.cols <- pca.x.norm # tests prcomp way
eigenfaces.as.cols <- principal.component.scores # tests do it yourself way

#### Then run the rest

# goal image - plotted later:
to.reconstruct.file <- paste(WD, "/CroppedYale/yaleB05/yaleB05_P00A+010E+00.pgm", sep="")
to.reconstruct.v <- as.vector(getChannels(read.pnm(to.reconstruct.file)))

top.eigenfaces <- lapply(c(1:24), function(i) { ith.eigenface.for.image.v(i, to.reconstruct.v, eigenfaces.as.cols) })
top.eigenfaces.cum <- cumsum.eigenfaces(top.eigenfaces)

# finally, plotting the evidence
par(mfrow = c(5,5))
for (i in seq_along(top.eigenfaces.cum)) {
  plot.row.v(top.eigenfaces.cum[[i]], paste("Sum of", i, "eigenfaces"))
}
plot.row.v(to.reconstruct.v, "Target Face 1e")