Singular Value Decomposition and Image Processing

# cleanup
rm(list=ls())

dat <- seq(1,36,2)

Y <- matrix(dat,ncol=6)
Y

##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    1    7   13   19   25   31
## [2,]    3    9   15   21   27   33
## [3,]    5   11   17   23   29   35

# Apply SVD to get U, V, and D
s <- svd(Y)
U <- s$u
V <- s$v
D <- diag(s$d) ##turn it into a matrix

Yhat <- U %*% D %*% t(V)
Yhat

##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    1    7   13   19   25   31
## [2,]    3    9   15   21   27   33
## [3,]    5   11   17   23   29   35

resid <- Y - Yhat
max(abs(resid))

## [1] 2.309264e-14

# source("http://bioconductor.org/biocLite.R")
# biocLite("ripa", dependencies=TRUE)
# biocLite("rARPACK", dependencies=TRUE)
# install.packages('devtools')
# library(devtools)
# install_github('ririzarr/rafalib')

# Load libraries
library(rARPACK)
library(ripa)    #  function "imagematrix"
library(EBImage)
library(jpeg)
library(png)
#
library(rafalib)
mypar2()
#
# Read the image
img <- readImage("images/pansy.jpg") 
dim(img)

## [1] 465 600   3

display(img, method = "raster")

# Convert image to greyscale
r <- imagematrix(img, type = "grey")

# display(r, method = "raster", all=TRUE)
display(r, method = "raster")

dim(r)

## [1] 465 600

str(r)

##  imagematrix [1:465, 1:600] 0.345 0.282 0.29 0.306 0.29 ...
##  - attr(*, "type")= chr "grey"
##  - attr(*, "class")= chr [1:2] "imagematrix" "matrix"

# Apply SVD to get u, v, and d
r.svd <- svd(r)

u <- r.svd$u
v <- r.svd$v
d <- diag(r.svd$d)
dim(d)

## [1] 465 465

## [1] 465 465

# check svd$d values 
# Plot the magnitude of the singular values
sigmas = r.svd$d # diagonal matrix (the entries of which are known as singular values)
plot(1:length(r.svd$d), r.svd$d, xlab="i-th r.svd$d", ylab="r.svd$d",  main="Singular Values");

plot(1:length(r.svd$d), cumsum(r.svd$d) / sum(r.svd$d), main="Cumulative Percent of Total Sigmas");

# first approximation
u1 <- as.matrix(u[-1, 1])
v1 <- as.matrix(v[-1, 1])
d1 <- as.matrix(d[1, 1])
l1 <- u1 %*% d1 %*% t(v1)
l1g <- imagematrix(l1, type = "grey")
#plot(l1g, useRaster = TRUE)
display(l1g, method = "raster", all=TRUE)

# more approximation
depth <- 5
us <- as.matrix(u[, 1:depth])
vs <- as.matrix(v[, 1:depth])
ds <- as.matrix(d[1:depth, 1:depth])
ls <- us %*% ds %*% t(vs)
lsg <- imagematrix(ls, type = "grey")

## Warning: Pixel values were automatically clipped because of range over.

#plot(lsg, useRaster = TRUE)
display(lsg, method = "raster")

# more approximation
depth <- 20
us <- as.matrix(u[, 1:depth])
vs <- as.matrix(v[, 1:depth])
ds <- as.matrix(d[1:depth, 1:depth])
ls <- us %*% ds %*% t(vs)
lsg <- imagematrix(ls, type = "grey")

## Warning: Pixel values were automatically clipped because of range over.

#plot(lsg, useRaster = TRUE)
display(lsg, method = "raster")

factorize = function(m, k){
  r = svds(m[, , 1], k);
  g = svds(m[, , 2], k);
  b = svds(m[, , 3], k);
  return(list(r = r, g = g, b = b));
}

recoverimg = function(lst, k){
  recover0 = function(fac, k){
    dmat = diag(k);
    diag(dmat) = fac$d[1:k];
    m = fac$u[, 1:k] %*% dmat %*% t(fac$v[, 1:k]);
    m[m < 0] = 0;
    m[m > 1] = 1;
    return(m);
  }

  r = recover0(lst$r, k);
  g = recover0(lst$g, k);
  b = recover0(lst$b, k);
  m = array(0, c(nrow(r), ncol(r), 3));
  m[, , 1] = r;
  m[, , 2] = g;
  m[, , 3] = b;
  return(m);
}

rawimg <- readJPEG("images/pansy.jpg");

lst = factorize(rawimg, 100);

neig = c(1, 5, 20, 50, 100);
for(i in neig){
  m = recoverimg(lst, i);
  writeJPEG(m, sprintf("images/svd_%d.jpg", i), 0.95);
  #display(m, method = "raster")
  fname <- sprintf("images/svd_%d.jpg", i)
  display(readImage(fname), title="svd_%d", method = "raster")
}