Maciej Czajka 85ba70d8b3 add calculate function

2022-05-18 00:08:40 +02:00

24 KiB

Raw Blame History

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
import scipy.linalg as la
from PIL import Image
from ipywidgets import interact

# zmień szerokość komórki
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:80% !important; }</style>"))

/var/folders/tq/jq5nwbnj7v10tls99x99qbh40000gn/T/ipykernel_57719/1134982733.py:12: 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

#Ta linia jest wymagana do wyświetlania wizualizacji w przeglądarce
%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_4u_4v_4^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 calculate(A):
    m = A.shape[0]
    n = A.shape[1]
    S = np.zeros(n)

    # finding eigenvectors with biggest eigenvalues of A*transpose(A)
    helper = np.matmul(A, np.transpose(A))
    eigenvalues, eigenvectors = la.eigh(helper)
    # descending sort of all the eigenvectors according to their eigenvalues
    index = eigenvalues.argsort()[::-1]
    eigenvalues = eigenvalues[index]
    eigenvectors = eigenvectors[:, index]
    U = eigenvectors

    # S is a diagonal matrix that keeps square root of eigenvalues
    j = 0
    for i in eigenvalues:
        if j == S.size:
            break
        elif i >= 0:
            S[j] = np.sqrt(i)
            j += 1
    # same finding process for transpose(A)*A
    helper = np.matmul(np.transpose(A), A)
    eigenvalues, eigenvectors = la.eigh(helper)
    # descending sort of all the eigenvectors according to their eigenvalues
    index = eigenvalues.argsort()[::-1]
    eigenvalues = eigenvalues[index]
    eigenvectors = eigenvectors[:, index]
    V = np.transpose(eigenvectors)

    # sorting S in descending order
    S[::-1].sort()
    # print_to_file(S)

    return U, S, V

def compress_svd(image,k):
    """
    Wykonaj dekompozycję SVD, a następnie okrojoną rekonstrukcję (przy użyciu k wartości/wektorów osobliwych)
    """
#     U,s,V = svd(image,full_matrices=False)
    U,s,V = calculate(image)
    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):
    """
    Kompresuje obrazy w skali szarości i wyświetla zrekonstruowany obraz.
    Wyświetla również wykres wartości singularnych
    """
    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(100 - compression_ratio)+"%")
    axes[1].imshow(reconst_img,cmap='gray')
    axes[1].axis('off')
    fig.tight_layout()
    # compression rate = 100% * (k * (height + width + k)) / (height + width)

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):
  """
    Funkcja do obliczania maksymalnej wartości zakresu suwaka "k"
  """
  img = gray_images[img_name]
  m,n = img.shape
  return m*n/(m+n+1)

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

def print_matrices(img_name, k):
    """
    Wyświetlanie macierzy U V S wraz z wymiarami.
    """

    if img_name not in gray_images.keys():
        return
    
    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)
    U,s,V = calculate(image)
    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', 'coffee', 'rocket', 'koala', '…

<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')),
    "teacher": img_as_float(Image.open('teacher.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):
    """
    Kompresowanie i wyświetlanie zrekonstruowanego obrazu kolorowego przy użyciu metody reshape.
    """
    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(100 - 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')

def print_color_matrices(img_name,k):
    """
    Wyświetlanie macierzy U V S wraz z wymiarami.
    """
    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)

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

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

<function __main__.print_color_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', 'rocket', 'koala', '…

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):
    """
    Kompresowanie i wyświetlanie zrekonstruowanego obrazu kolorowego przy użyciu metody warstwowej.
    """
    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(100- 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(print_color_matrices,img_name=list_widget,k=int_slider_widget);

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

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

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

24 KiB Raw Blame History