26 KiB
import pandas as pd
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
from numpy.linalg import eig
from math import isclose
import pprint
# zmień szerokość komórki
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:80% !important; }</style>"))
/var/folders/8w/3c34c7kd2n144tvm764_pdbw0000gq/T/ipykernel_3019/1657712203.py:16: 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 inlineKompresja 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żywamyNumPy 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ścireconst_matrix = np.dot(U[:,:k],np.dot(np.diag(s[:k]),V[:k,:]))
def calculate(A):
n = A.shape[0]
m = A.shape[1]
S = np.zeros(n)
helper = np.matmul(A, np.transpose(A))
eigenvalues, eigenvectors = eig(helper)
index = eigenvalues.argsort()[::-1]
eigenvalues = np.real(eigenvalues)
eigenvectors = np.real(eigenvectors)
eigenvalues = eigenvalues[index]
eigenvectors = eigenvectors[:, index]
U = eigenvectors
helper2 = np.matmul(np.transpose(A), A)
eigenvalues2, eigenvectors2 = eig(helper2)
index2 = eigenvalues2.argsort()[::-1]
eigenvalues2 = np.real(eigenvalues2)
eigenvectors2 = np.real(eigenvectors2)
eigenvalues2 = eigenvalues2[index2]
eigenvectors2 = eigenvectors2[:, index2]
V = np.transpose(eigenvectors2)
j = 0
for i in eigenvalues2:
if j == S.size:
break
elif i >= 0:
S[j] = np.sqrt(i)
j += 1
S[::-1].sort()
return U, S, Vdef check_sign_V(V, builtin):
if builtin.shape[0] < builtin.shape[1]:
for i in range(0,builtin.shape[0]):
for j in range(0,builtin.shape[1]):
if builtin[j][i] < 0.0 and V[j][i] > 0.0:
V[j][i] *= -1
elif builtin[j][i] > 0.0 and V[j][i] < 0.0:
V[j][i] *= -1
else:
for i in range(0,builtin.shape[0]):
for j in range(0,builtin.shape[1]):
if builtin[i][j] < 0.0 and V[i][j] > 0.0:
V[i][j] *= -1
elif builtin[i][j] > 0.0 and V[i][j] < 0.0:
V[i][j] *= -1
return V
def check_sign_U(U, builtin):
if builtin.shape[0] < builtin.shape[1]:
for i in range(0,builtin.shape[0]):
for j in range(0,builtin.shape[1]):
if builtin[j][i] < 0.0 and U[j][i] > 0.0:
U[j][i] *= -1
elif builtin[j][i] > 0.0 and U[j][i] < 0.0:
U[j][i] *= -1
else:
for i in range(0,builtin.shape[0]):
for j in range(0,builtin.shape[1]):
if builtin[i][j] < 0.0 and U[i][j] > 0.0:
U[i][j] *= -1
elif builtin[j][j] > 0.0 and U[i][j] < 0.0:
U[i][j] *= -1
return Udef compress_svd(image,k):
"""
Wykonaj dekompozycję SVD, a następnie okrojoną rekonstrukcję (przy użyciu k wartości/wektorów osobliwych)
"""
image = np.real(image)
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,sPrzykładowa macierz
a = np.array([[1, 0, 0, 0, 2],
[0, 0, 3, 0, 0],
[0, 0, 0, 0, 0],
[0, 4, 0, 0, 0]])
U,s,V = calculate(a)
U1,s1,V1 = svd(a)
U = check_sign_U(U, U1)
V = check_sign_V(V, V1)
U_dataframe = pd.DataFrame(U)
s_dataframe = pd.DataFrame(s)
V_dataframe = pd.DataFrame(V)
print("Matrix U:")
print(U_dataframe)
print("\nMatrix V:")
print(V_dataframe)
print("\nMatrix s:")
print(s_dataframe)
print('\n--------------------------------------\n')
k = 4
reconst_matrix4 = np.dot(U[:,:k],np.dot(np.diag(s[:k]),V[:k,:]))
recon_dataframe_matrix = pd.DataFrame(reconst_matrix4)
print('Reconstructed matrix: \n')
print(recon_dataframe_matrix)Matrix U:
0 1 2 3
0 0.0 0.0 1.0 0.0
1 0.0 1.0 0.0 0.0
2 0.0 0.0 0.0 -1.0
3 1.0 0.0 0.0 0.0
Matrix V:
0 1 2 3 4
0 0.000000 1.0 0.0 0.0 0.000000
1 0.000000 0.0 1.0 0.0 0.000000
2 0.447214 0.0 0.0 0.0 0.894427
3 0.000000 0.0 0.0 1.0 0.000000
4 -0.894427 0.0 0.0 0.0 0.447214
Matrix s:
0
0 4.000000
1 3.000000
2 2.236068
3 0.000000
--------------------------------------
Reconstructed matrix:
0 1 2 3 4
0 1.0 0.0 0.0 0.0 2.0
1 0.0 0.0 3.0 0.0 0.0
2 0.0 0.0 0.0 0.0 0.0
3 0.0 4.0 0.0 0.0 0.0
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)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')W celu zbadania, jak jakość zrekonstruowanego obrazu zmienia się wraz z $k$ należy użyć poniższego interaktywnego widżetu.
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]}\n")
U,s,V = svd(image,full_matrices=False)
U_dataframe = pd.DataFrame(U[:,:k])
print(f"U MATRIX:\n {U_dataframe}\n")
print('\n', '*' * 100, '\n')
print(f"\nShape of S matrix: {s[:k].shape[0]}\n")
s_dataframe = pd.DataFrame(np.diag(s[:k]))
print(f"S MATRIX:\n {s_dataframe}")
print('\n', '*' * 100, '\n')
V_dataframe = pd.DataFrame(V[:k,:].T)
print(f"V MATRIX:\n {V_dataframe}")
interact(print_matrices, img_name=list_widget, k=int_slider_widget)interactive(children=(Dropdown(description='img_name', options=('cat', 'astro', 'line', 'camera', 'coin', 'clo…<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', 'line', 'camera', 'coin', 'clo…Ł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)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_color = 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_color.value
int_slider_widget.max = compute_k_max_color_images(img_name)
list_widget_color.observe(update_k_max_color,'value')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.
interact(compress_show_color_images_reshape,img_name=list_widget_color,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)]
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_color = 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_color.value
int_slider_widget.max = compute_k_max_color_images_layers(img_name)
list_widget_color.observe(update_k_max_color_layers,'value')interact(compress_show_color_images_layer,img_name=list_widget_color,k=int_slider_widget);interactive(children=(Dropdown(description='img_name', options=('cat', 'astro', 'coffee', 'rocket', 'koala', '…