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
from ipywidgets import interact
# change the cell width
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:80% !important; }</style>"))

/var/folders/t3/dwnz0lh916ng4w7bzf0z56ym0000gn/T/ipykernel_20610/17056051.py:11: DeprecationWarning: Importing display from IPython.core.display is deprecated since IPython 7.14, please import from IPython display
  from IPython.core.display import display, HTML

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

Kompresja obrazów z wykorzystaniem rozkładu SVD

Metodę SVD zastosujemy do kompresji obrazów. SVD dokonuje dekompozycji macierzy prostokątnej $M$ na trzy części. $M=U\Sigma V^T$ -

• $U$ - macierz, w której kolumny są wektorami własnymi macierzy $MM^T$ (lewe wektory osobliwe - 'left singular vectors'). Wektory te tworzą układ ortonormalny.
• $\Sigma$ - macierz diagonalna, która na swojej diagonalii ma nieujemne wartości osobliwe (szczególne), czyli pierwiastki wartości własnych macierzy $M^TM$ uporządkowane nierosnąco
• $V$ - macierz, w której kolumny są wektorami własnymi macierzy $M^TM$ (prawe wektory osobliwe - 'right singular vectors'). Wektory te tworzą układ ortonormalny.

SVD polega na rekonstrukcji oryginalnej macierzy jako kombinacji liniowej kilku macierzy rangi jeden. Macierz rangi jeden można wyrazić jako iloczyn zewnętrzny dwóch wektorów kolumnowych.

$M=\sigma_1u_1v_1^T+\sigma_2u_2v_2^T+\sigma_3u_3v_3^T+\sigma_3u_3v_3^T+....$ .

$rank M=r$ . $M=\sum_{i=1}^{r} \sigma_iu_iv_i^T$

Macierz o randze r będzie miała r takich wyrazów.

Tutaj $\sigma_1,\sigma_2,\sigma_3 ...$ są wartościami osobliwymi. $u_1,u_2,u_3 ...$ i $v_1,v_2,v_3 ...$ są kolejnymi kolumnami (wektorami własnymi) z macierzy U i V.

Kompresja obrazu przy użyciu SVD polega na wykorzystaniu faktu, że tylko nieliczne wartości osobliwe są duże. Chociaż obrazy ze świata rzeczywistego mają pełną rangę, to ich efektywna ranga jest niska, co oznacza, że tylko kilka wartości osobliwych rozkładu SVD obrazów będzie dużych.

skimage - biblioteka do przetwarzania obrazów

Do pracy z obrazami w Pythonie używamy biblioteki do przetwarzania obrazów skimage (z rodziny pakietów sci-kit). Ma ona moduł o nazwie data, który udostępnia zestaw obrazów. Wczytujemy kilka obrazów i konwertujemy je do formatu skali szarości. Obrazy te są przechowywane w słowniku 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(),
        "text":data.text(),
        "page":data.page(),
        "blobs":data.binary_blobs(),
        "coffee":rgb2gray(img_as_float(data.coffee()))
}

SVD w python'ie

Używamy funkcji svd z biblioteki numpy.linalg do obliczenia rozkładu SVD naszej macierzy. Funkcja svd zwraca U,s,V .

• U macierz, w której kolumny są wektorami własnymi macierzy $MM^T$

• s jest tablicą numpy rangi 1 z wartościami osobliwymi

• V macierz, w której wiersze są wektorami własnymi macierzy $M^TM$ - odpowiednik $V^T$ w tradycyjnej literaturze algebry liniowej

  Zrekonstruowana aproksymacja oryginalnej macierzy jest wykonywana przy użyciu podzbioru wektorów osobliwych, jak poniżej w funkcji compress_svd. Używamy NumPy Array Slicing, aby wybrać k wektorów osobliwych i wartości osobliwych. Zamiast przechowywać $m\times n$ wartości dla oryginalnego obrazu, możemy teraz przechowywać $k(m+n)+k$ wartości

    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

Kompresja obrazów w skali szarości

Poniższa funkcja compress_show_gray_images przyjmuje nazwę obrazu (img_name) i liczbę wartości/wektorów osobliwych (k), które mają być użyte w skompresowanej rekonstrukcji. Na wykresie są wartości osobliwe oraz obraz.

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()

W celu zbadania, jak jakość zrekonstruowanego obrazu zmienia się wraz z $k$ należy użyć poniższego interaktywnego widżetu.

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)

interactive(children=(Dropdown(description='img_name', options=('cat', 'astro', 'camera', 'coin', 'clock', 'te…

<function __main__.print_matrices(img_name, k)>

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

interactive(children=(Dropdown(description='img_name', options=('cat', 'astro', 'camera', 'coin', 'clock', 'te…

Ładowanie kolorowych obrazów

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

Kompresja kolorowych obrazów

Obrazy kolorowe są reprezentowane w Pythonie jako trójwymiarowe tablice numpy - trzeci wymiar reprezentuje wartości kolorów (czerwony, zielony, niebieski). Jednak metoda svd ma zastosowanie do macierzy dwuwymiarowych. W tym celu musimy przekonwertować tablicę trójwymiarową na tablicę dwuwymiarową, zastosować svd i zrekonstruować ją z powrotem jako tablicę trójwymiarową. Istnieją dwa sposoby, aby to zrobić:

• Metoda reshape
• Metoda layers (warstwowa)

Obie te metody pokażemy poniżej.

Metoda reshape do kompresji kolorowych obrazów

Metoda ta polega na spłaszczeniu trzeciego wymiaru tablicy obrazów do drugiego wymiaru przy użyciu metody reshape z biblioteki numpy.

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

Dekompozycja svd jest stosowana do wynikowej, przekształconej tablicy i rekonstruowana z żądaną liczbą wartości/wektorów osobliwych. Tablica obrazów jest przekształcana z powrotem do trzech wymiarów przez kolejne wywołanie metody reshape.

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)

Oto interaktywny widżet do badania kompresji obrazów kolorowych metodą reshape. Przeciągając suwak w celu zmiany wartości $k$, można zaobserwować, jak zmienia się jakość obrazu. Można także badać różne obrazy, wybierając je za pomocą rozwijanego widżetu.

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)

interactive(children=(Dropdown(description='img_name', options=('cat', 'astro', 'coffee', 'koala', 'orange'), …

<function __main__.print_matrices(img_name, k)>

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

interactive(children=(Dropdown(description='img_name', options=('cat', 'astro', 'coffee', 'koala', 'orange'), …

Metoda layers (warstwowa) do kompresji kolorowych obrazów

W funkcji compress_show_color_images_layer, traktujemy kolorowy obraz jako stos 3 oddzielnych dwuwymiarowych obrazów (warstwy czerwona, niebieska i zielona). Stosujemy okrojoną rekonstrukcję svd na każdej dwuwymiarowej warstwie osobno.

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

I ponownie łączymy zrekonstruowane warstwy razem.

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)

Oto widżet do badania metody layers (warstwowej) kompresji obrazów kolorowych.

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);

interactive(children=(Dropdown(description='img_name', options=('cat', 'astro', 'coffee', 'koala', 'orange'), …