Tuesday, July 11, 2017

Image Compression Using SVD with R




I just finished the great Linear Algebra course from Prof. Gilbert Strang thanks to MIT OpenCourseWare. It was a fantastic series of lectures. The reason I started to that course was understanding Singular Value Decomposition at the first place. But I learned a lot throughout the course. I definitely recommend you to watch some of them even if you have already taken a course on Linear Algebra. After finishing the course, I decided to write a blog using concepts that I learned from that course. Let's start!


Singular Value Decomposition is another matrix factorization. After performing SVD  to the matrix we get 3 different matrices that are orthogonal, diagonal, orthogonal. Let transpose of Q equal to Q'. Eigenvalue - eigenvector factorization QΛQ' is an example of SVD. It has orthonormal, diagonal and orthonormal matrices(we can convert orthogonal to orthonormal by simply dividing with the length of the vector). But not all matrices can be factorized into that great form. If matrix A is not symmetric positive definite we cannot get the eigenvalue - eigenvector decomposition. We need to look at  AA' and A'A and then get their eigenvectors.

So with matrix A, what we are looking for is an orthonormal basis. SVD factors the matrix into three matrices A =UDV'. U is orthonormal basis in the column space and V is orthonormal basis in the row space. Matrix U is eigenvectors of AA' and matrix V is eigenvectors of A'A. Eigenvalues of AA' and A'A are equal.

After applying those steps we can get our SVD factorization. But why does it necessary? Where we can use it? In SVD factorization, we are getting a diagonal matrix D that has square roots of eigenvalues of A. The matrices that has the eigenvalues itself are AA' and A'A. So, when it comes to matrix A, those square roots of eigenvalues are our singular values. And those singular values are varying. Some of them are significant and others are extremely small. If we just get some significant singular values we can explain or show our matrix to some extend. This concept can be applied to image compression. If we have an image that contains 512 x 512 pixels we can compress it into less pixels using SVD. We will do that by selecting low level ranked matrices.

Let's apply this using R.
JPEG package can be used to read and write jpg files. I will use a bacteria image to compress.




> library(jpeg)
> bacteria <- readJPEG('bac.jpg')
> ncol(bacteria)
> [1] 600
> nrow(bacteria)
> [1] 431
> dim(bacteria)
> [1] 431 600   3
RGB color model uses red, green and blue colors. As can be seen from the dimensions, our image has 431x600 pixels for each of the color types. We can split them into three different matrixes.
> r <- bacteria[,,1] 
> g <- bacteria[,,2]
> b <- bacteria[,,3] 
Now we have our RGB colors. We are ready to apply SVD.
> bacteria.r <- svd(r)
> bacteria.g <- svd(g)
> bacteria.b <- svd(b)
Now lets make a list using them.
> rgb.list <- list(bacteria.r, bacteria.g, bacteria.b)
In that rgb.list, we have 3 lists and they each also another list of 3. Each of them has singular value decomposition results inside. We need to reach those values inside that list. For this reason, sapply function that allows to reach variables features inside the list is useful.
> comp <- sapply(rgb.list, function(i){
> compressed = i$u[,1:3] %*% diag(i$d[1:3]) %*% t(i$v[,1:3])
  }, simplify = 'array')
Here we select only the first 3 components of SVD instead of getting all 431. In other words, we have a rank 3 matrix. Full rank would be our original image.
> writeJPEG(comp, paste('compressed/','bacteria_svd_level_', 3, '.jpg' ))
Here is our rank 3




This is the 53 rank version of the same processes. Only code that changed is [,1:3] to [,1:53].





With 53 rank we can get the most of the image. Even less can be enough to see those bacterias. SVD helped us to decrease the size of the image. It is helpful when we are working with lots of images. It makes it easier to store and retrieve the data.









No comments:

Post a Comment

Gibbs Sampler

Gibbs Sampler In this blog we are going to discuss Gibbs Sampler with an applied example. But we will not dive into Gibbs Sampler direc...