"A notebook to demonstrate the task of image compression using singular value decomposition. The core idea is to express the image(a matrix) as a product of other matrices which can be computed or derived from the source matrix. This then gives us the benefit to analyze how much of the matrix is worth keeping for the sake of a near-perfect reconstruction."
from google.colab import drive
drive.mount('/gdrive')
%cd /gdrive
Mounted at /gdrive /gdrive
from matplotlib.image import imread
import matplotlib.pyplot as plt
import numpy as np
import os
import matplotlib
matplotlib.rcParams['figure.facecolor'] = '#ffffff'
plt.rcParams['figure.figsize'] = [14,14]
The meaning of gray = np.mean(rgb, -1) is the same as:
r, g, b = rgb[:,:,0], rgb[:,:,1], rgb[:,:,2]
gray = 0.2989 * r + 0.5870 * g + 0.1140 * b
Note: The image used here is one of my own designs, please do not use it freely.
I =imread('/gdrive/My Drive/Colab Notebooks/images/byomkesh.jpg')
I_gray=np.mean(I,-1)
img =plt.imshow(I_gray)
img.set_cmap('gray')
plt.axis('off')
plt.show()
Singular value decomposition takes a rectangular matrix ($A_{n*p}$), and states the following:
$A_{n*p} = U_{n*n}*S_{n*p}*V^{T}_{p*p}$
where $U$ and $V$ are orthogonal.
Calculating the SVD consists of finding the eigenvalues and eigenvectors of $AA^T$ and $A^TA$. The eigenvectors of $A^TA$ make up the columns of $V$ , the eigenvectors of $AA^T$ make up the columns of $U$. Also, the singular values in $S$ are square roots of eigenvalues from $AA^T$ or $A^TA$. The singular values are the diagonal entries of the S matrix and are arranged in descending order.
The singular values are always real numbers. If the matrix $A$ is a real matrix, then $U$ and $V$ are also real.
We use the np.linalg.svd() function here. The images are reconstructed at different ranks from 5 to 2500 to visualise how much of the compression works
If we look closely at the code below we will see that we use the first r columns of U and the first r rows of VT and the rXr blocks of S to reconstruct the original matrix or the image. Now this is what is known as a rank-approximation where we iteratively reconstruct the image using different values of r.
#using the np.linalg.svd() function to compute SVD
U, S, VT = np.linalg.svd(I_gray,full_matrices=False)
S=np.diag(S)
#the following for loop will calculate
#the reconstruction of the original image using
#a total of 6 ranks ranging from 5 to 2500
j=0
for r in (5, 20, 100, 500, 1500, 2500):
#Construct approx images
Iapprox = U[:,:r]@S[0:r,:r]@VT[:r,:]
plt.subplot(2,3,j+1)
j+=1
img =plt.imshow(Iapprox)
img.set_cmap('gray')
plt.axis('off')
plt.title('r='+str(r))
plt.show()
This is done to visualise on a log scale how much of the matrix is useful or relevant to reconstruction and how much is redundant. So we can see that there is almost a flat line after 1500. So discarding the values from 1500 to the end will still give us a good enough reconstruction.
#plotting the singular values ina log scale
plt.figure(1)
plt.semilogy(np.diag(S))
plt.title('Singular Values')
plt.show()
#plotting the cumulative sum of the singular values
plt.figure(2)
plt.plot(np.cumsum(np.diag(S))/np.sum(np.diag(S)))
plt.title('Singular Values Cumulative Sum')
plt.show()
#final compression of the image
rnew = 1500
Ifinal = U[:,:rnew]@S[0:rnew,:rnew]@VT[:rnew,:]
img =plt.imshow(Ifinal)
img.set_cmap('gray')
plt.axis('off')
plt.show()
This notebook and subsequent blog post has been significantly inspired from the following resources: