svd_mpsic/jupyter.ipynb

import IPython
import numpy as np
import matplotlib.pyplot as plt
from skimage import data
from skimage.color import rgb2gray
from skimage import img_as_ubyte,img_as_float
from numpy.linalg import svd
from PIL import Image
# change the cell width
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:80% !important; }</style>"))

#This line is required to display visualizations in the browser
%matplotlib inline

SVD and Image compression

Now we will explore how to apply Singular Value Decomposition of a matrix to the problem of image compression. SVD decomposes a rectangular matrix $M$ to a three parts. $M=U\Sigma V^T$ -

• $U$ - matrix of left singular vectors in the columns
• $\Sigma$ - diagonal matrix with singular values
• $V$ - matrix of right singular vectors in the columns

SVD in effect involves reconstructing the original matrix as a linear combination of several rank one matrices. A rank one matrix can be expressed as a outer product of two column vectors.

$M=\sigma_1u_1v_1^T+\sigma_2u_2v_2^T+\sigma_3u_3v_3^T+\sigma_3u_3v_3^T+....$ . A matrix of rank r will have r terms of these.

Here $\sigma_1,\sigma_2,\sigma_3 ...$ are singular values. $u_1,u_2,u_3 ...$ and $v_1,v_2,v_3 ...$ are left and right singular vectors respectively.

Image compression using SVD involves taking advantage of the fact that very few of the singular values are large. Although images from the real world are of full rank, they have low effective rank which means that only few of the singular values of the SVD of images will be large.

skimage image processing library

We will use skimage image processing library (from sci-kit family of packages) for working with images in python. skimage has a module called data which makes available a set of images for exploration. We will load some images and convert them into a gray scale format. These images are stored in a python dict object gray_images.

gray_images = {
        "cat":rgb2gray(img_as_float(data.chelsea())),
        "astro":rgb2gray(img_as_float(data.astronaut())),
        "camera":data.camera(),
        "coin": data.coins(),
        "clock":data.clock(),
        "blobs":data.binary_blobs(),
        "coffee":rgb2gray(img_as_float(data.coffee()))
}

svd in python

We will use numpy.linalg library's svd function to compute svd of a matrix in python. The svd function returns U,s,V .

• U has left singular vectors in the columns

• s is rank 1 numpy array with singular values

• V has right singular vectors in the rows -equivalent to $V^T$ in traditional linear algebra literature

  The reconstructed approximation of the original matrix is done using a subset of singular vectors as below in the compress_svd function . We use numpy array slicing to select k singular vectors and values. Instead of storing $m\times n$ values for the original image, we can now store $k(m+n)+k$ values

    reconst_matrix = np.dot(U[:,:k],np.dot(np.diag(s[:k]),V[:k,:]))
    
  

def compress_svd(image,k):
    """
    Perform svd decomposition and truncated (using k singular values/vectors) reconstruction
    returns
    --------
      reconstructed matrix reconst_matrix, array of singular values s
    """
    U,s,V = svd(image,full_matrices=False)
    reconst_matrix = np.dot(U[:,:k],np.dot(np.diag(s[:k]),V[:k,:]))

    return reconst_matrix,s

Compress gray scale images

The function compress_show_gray_images below takes in the image name (img_name) and number of singular values/vectors(k) to be used in the compressed reconstruction. It also plots the singular values and the image.

def compress_show_gray_images(img_name,k):
    """
     compresses gray scale images and display the reconstructed image.
     Also displays a plot of singular values
    """
    image=gray_images[img_name]
    original_shape = image.shape
    reconst_img,s = compress_svd(image,k)
    fig,axes = plt.subplots(1,2,figsize=(8,5))
    axes[0].plot(s)
    compression_ratio =100.0* (k*(original_shape[0] + original_shape[1])+k)/(original_shape[0]*original_shape[1])
    axes[1].set_title("compression ratio={:.2f}".format(compression_ratio)+"%")
    axes[1].imshow(reconst_img,cmap='gray')
    axes[1].axis('off')
    fig.tight_layout()

Use the below interactive widget to explore how the quality of the reconstructed image varies with $k$

def compute_k_max(img_name):
  """
    utility function for calculating max value of the slider range
  """
  img = gray_images[img_name]
  m,n = img.shape
  return m*n/(m+n+1)

#set up the widgets
import ipywidgets as widgets

list_widget = widgets.Dropdown(options=list(gray_images.keys()))
int_slider_widget = widgets.IntSlider(min=1,max=compute_k_max('cat'))
def update_k_max(*args):
  img_name=list_widget.value
  int_slider_widget.max = compute_k_max(img_name)
list_widget.observe(update_k_max,'value')

# Print matrices
def print_matrices(img_name, k):
    image=gray_images[img_name]
    original_shape = image.shape
    print(f"Input image dimensions. Width:{original_shape[1]} Height:{original_shape[0]}")

    U,s,V = svd(image,full_matrices=False)
    print(f"Shape of U matrix: {U[:,:k].shape}")
    print(f"U MATRIX: {U[:,:k]}")
    print('*' * 100)
    print(f"Shape of S matrix: {s[:k].shape}")
    print(f"S MATRIX: {np.diag(s[:k])}")
    print('*' * 100)
    print(f"Shape of V matrix: {V[:k,:].shape}")
    print(f"V MATRIX: {V[:k,:]}")

interact(print_matrices, img_name=list_widget, k=int_slider_widget)

interact(compress_show_gray_images,img_name=list_widget,k=int_slider_widget);

Load color images

color_images = {
    "cat":img_as_float(data.chelsea()),
    "astro":img_as_float(data.astronaut()),
    "coffee":img_as_float(data.coffee()),
    "koala": img_as_float(Image.open('koala.jpeg')),
    "orange": img_as_float(Image.open('orange.jpeg'))
}

Compress color images

Color images are represented in python as 3 dimensional numpy arrays - the third dimension to represent the color values (red,green blue). However, svd method is applicable to two dimensional matrices. So we have to find a way to convert the 3 dimensional array to 2 dimensional arrays, apply svd and reconstruct it back as a 3 dimensional array . There are two ways to do it. We will show both these methods below .

• reshape method
• Layer method

Reshape method to compress a color image

This method involves flattening the third dimension of the image array into the second dimension using numpy's reshape method .

image_reshaped = image.reshape((original_shape[0],original_shape[1]*3))

The svd decomposition is applied on the resulting reshaped array and reconstructed with the desired number of singular values/vectors. The image array is reshaped back to the three dimensions by another call to reshape method.

image_reconst = image_reconst.reshape(original_shape)

def compress_show_color_images_reshape(img_name,k):
    """
     compress and display the reconstructed color image using the reshape method 
    """
    image = color_images[img_name]
    original_shape = image.shape
    image_reshaped = image.reshape((original_shape[0],original_shape[1]*3))
    image_reconst,_ = compress_svd(image_reshaped,k)
    image_reconst = image_reconst.reshape(original_shape)
    compression_ratio =100.0* (k*(original_shape[0] + 3*original_shape[1])+k)/(original_shape[0]*original_shape[1]*original_shape[2])
    plt.title("compression ratio={:.2f}".format(compression_ratio)+"%")
    plt.imshow(image_reconst)

Here is the interactive widget to explore image compression of color images using the reshape method. By dragging the slider to vary $k$, observe how image quality varies. Also, we can explore different images by selecting through the drop down widget.

def compute_k_max_color_images(img_name):
  image = color_images[img_name]
  original_shape = image.shape
  return (original_shape[0]*original_shape[1]*original_shape[2])//(original_shape[0] + 3*original_shape[1] + 1)


list_widget = widgets.Dropdown(options=list(color_images.keys()))
int_slider_widget = widgets.IntSlider(min=1,max=compute_k_max_color_images('cat'))
def update_k_max_color(*args):
  img_name=list_widget.value
  int_slider_widget.max = compute_k_max_color_images(img_name)
list_widget.observe(update_k_max_color,'value')

interact(print_matrices, img_name=list_widget, k=int_slider_widget)

interact(compress_show_color_images_reshape,img_name=list_widget,k=int_slider_widget);

Layers method to compress color images

In the function compress_show_color_images_layer, we treat a color image as a stack of 3 seperate two dimensional images (Red,blue and green layers) . We apply the truncated svd reconstruction on each two dimensional layer seperately.

image_reconst_layers = [compress_svd(image[:,:,i],k)[0] for i in range(3)]

And we put back the reconstructed layers together.

image_reconst = np.zeros(image.shape)
for i in range(3):
   image_reconst[:,:,i] = image_reconst_layers[i]

def compress_show_color_images_layer(img_name,k):
    """
     compress and display the reconstructed color image using the layer method 
    """
    image = color_images[img_name]
    original_shape = image.shape
    image_reconst_layers = [compress_svd(image[:,:,i],k)[0] for i in range(3)]
    print(image_reconst_layers)
    image_reconst = np.zeros(image.shape)
    for i in range(3):
        image_reconst[:,:,i] = image_reconst_layers[i]
    
    compression_ratio =100.0*3* (k*(original_shape[0] + original_shape[1])+k)/(original_shape[0]*original_shape[1]*original_shape[2])
    plt.title("compression ratio={:.2f}".format(compression_ratio)+"%")
    plt.imshow(image_reconst, vmin=0, vmax=255)

Here is the widget to explore layers method of compressing color images.

def compute_k_max_color_images_layers(img_name):
  image = color_images[img_name]
  original_shape = image.shape
  return (original_shape[0]*original_shape[1]*original_shape[2])// (3*(original_shape[0] + original_shape[1] + 1))


list_widget = widgets.Dropdown(options=list(color_images.keys()))
int_slider_widget = widgets.IntSlider(min=1,max=compute_k_max_color_images_layers('cat'))
def update_k_max_color_layers(*args):
  img_name=list_widget.value
  int_slider_widget.max = compute_k_max_color_images_layers(img_name)
list_widget.observe(update_k_max_color_layers,'value')

interact(compress_show_color_images_layer,img_name=list_widget,k=int_slider_widget);