Beeldcompressie met singulierewaardenontbinding

Engelse term: Singular Value Decomposition of SVD. Laad een plaatje als een Numpy-array A. Als het een kleuren-JPEG is, dan heeft de array drie dimensies: hoogte (rijen), breedte (kolommen) en kleurwaarden (R, G en B in een bepaalde volgorde die nu niet belangrijk is).

In [1]:
import numpy as np
import matplotlib.pyplot as plt

A = plt.imread("lion.jpg")
A.shape
Out[1]:
(2000, 3200, 3)

Bepaal de gemiddelde waarden van R, G, en B voor elke pixel, zodat het een grijswaardenplaatje wordt. De functie np.mean() neemt een matrix en middelt één van de dimensies ("axis") die beginnen te tellen bij nul. Hier dus de derde dimensie: axis=2 (of gebruik axis=-1 voor de laatste, wat in dit geval hetzelfde is).

In [2]:
B = np.mean(A, axis=2)
B.shape
Out[2]:
(2000, 3200)

Geef het plaatje weer met de berekende grijswaarden.

In [3]:
plt.figure(0, figsize=(20,20))
img = plt.imshow(B)
img.set_cmap("gray")
plt.axis("off")
plt.show()
2020-08-10T15:32:19.754549 image/svg+xml Matplotlib v3.3.0, https://matplotlib.org/

Bepaal de singulierewaardenontbinding met de functie np.linalg.svd(). Het resultaat bestaat uit twee matrices en de hoofddiagonaal van een matrix:

  • matrix U die de "eigenkolommen" van het beeld voorstelt,
  • matrix VT (V getransponeerd) die de "eigenrijen" van het beeld voorstelt,
  • hoofddiagonaal S1 (Sigma-ééndimensionaal) van de singulierewaardenmatrix Sigma.

Matrix S (Sigma) is een diagonale matrix, dat wil zeggen: elke waarde in S is nul behalve die op de hoofddiagonaal. De getallen op de diagonaal geven aan hoe belangrijk de kolommen in U en rijen in VT zijn. Het SVD-algoritme werkt zó, dat de waarden in S1 gesorteerd zijn van groot naar klein, wat betekent dat ook de eigenvectoren van het beeld (de kolommen in U en de rijen in VT) gesorteerd zijn op belangrijkheid.

Omdat de complete SVD enorme matrices genereert, gebruik je hier full_matrices=False om de grootte te reduceren, een methode die in Matlab de Economy SVD heet. Dit mag omdat Sigma in het algemeen niet vierkant is, terwijl alleen op de hoofddiagonaal waarden staan die niet nul zijn. Er is dus een rechthoekig stuk aan de rechter- of onderkant van Sigma met alleen maar nullen, dat je ongestraft weg kunt laten. Daardoor reduceert matrix U ook een stuk, wat anders alleen maar met nul vermenigvuldigd zou worden.

De originele beeldmatrix B is gegarandeerd exact gelijk aan het product van deze onderdelen van de SVD: B = U . S . VT, zelfs je de economy-versie neemt.

In [4]:
U, S1, VT = np.linalg.svd(B, full_matrices=False)  # calculate the "economy" SVD
S = np.diagflat(S1)                                # generate square diagonal matrix S from its main diagonal
[U.shape, S.shape, VT.shape]
Out[4]:
[(2000, 2000), (2000, 2000), (2000, 3200)]

Maak twee grafieken van de singuliere waarden, om te controleren of de ontbinding geschikt is om een gecomprimeerde versie van het beeld te maken.

Eerst de singuliere waarden zelf, op een logaritmische schaal. Als het goed is, is er een buigpunt: de waarden zakken eerst snel in, daarna langzamer. Dat betekent dat de eerste n eigenvectoren van het beeld zeer belangrijk zijn, en een groot gedeelte op het einde allemaal ongeveer even onbelangrijk. Daarna de genormaliseerde cumulatieve som van de waarden, waaraan je kunt zien hoeveel procent van alle beeldinformatie in de eerste n singuliere waarden zit.

Ter oriëntatie enkele verticale lijnen voor min of meer willekeurige waarden van n, en een horizontale lijn op 90%.

In [5]:
plt.figure(1, figsize=(20,7))

plt.subplot(121)
plt.vlines([50, 200, 500], S1[-1], S1[0], "r", "--")
plt.semilogy(S1)
plt.title("Singuliere waarden")

plt.subplot(122)
plt.hlines(0.9, 0, S1.size, "r", "--")
plt.vlines([50, 200, 500], 0, 1, "r", "--")
plt.plot(np.cumsum(S1) / np.sum(S1))
plt.title("Cumulatieve som")

plt.show()
2020-08-10T15:32:30.132898 image/svg+xml Matplotlib v3.3.0, https://matplotlib.org/

Nu waar het allemaal om draait: beeldcompressie door onbelangrijke informatie weg te gooien. Je reduceert de rang van de matrices U, S en VT door alleen de eerste n kolommen of rijen mee te nemen in de berekening van de benadering X:

  • reduceer het aantal kolommen van U tot n,
  • reduceer de vierkante diagonaalmatrix S tot n x n,
  • reduceer het aantal rijen van VT tot n.

Het komt goed uit dat de grootste waarden in S vooraan staan: daardoor bewaar je de belangrijke informatie en gooi je minder belangrijke informatie weg. Hoe belangrijk de waarden precies zijn, oftewel: bij welke n je kunt afkappen, analyseer je uit de singulierewaardengrafieken, of visueel uit het resultaat X van de benadering: hoe goed zijn onderstaande plaatjes?

In [6]:
plt.figure(2, figsize=(20,28))

i = 1
for n in (2, 5, 10, 20, 50, 100, 200, 500):
    X = U[:, :n] @ S[:n, :n] @ VT[:n, :]  # X approximates B by throwing away small singular values
    plt.subplot(4, 2, i)
    img = plt.imshow(X)
    img.set_cmap("gray")
    plt.axis("off")
    plt.title("n = " + str(n))
    i += 1

plt.show()
2020-08-10T15:32:33.242408 image/svg+xml Matplotlib v3.3.0, https://matplotlib.org/