SciPy: FFT, DCT a DST prováděná nad rastrovými obrázky

V předchozím článku o knihovně SciPy jsme se zabývali funkcemi určenými pro výpočet transformací jednorozměrných signálů. Jednalo se především o rychlou Fourierovu transformaci (FFT), což je algoritmus pro rychlý výpočet diskrétní Fourierovy transformace (DFT). Taktéž jsme se zmínili o diskrétní kosinové transformaci (DCT). Na tento článek dnes navážeme, protože si ukážeme zpracování 2D signálů, konkrétně rastrových obrázků. Budeme počítat rychlou dopřednou i zpětnou FFT, DCT i DST (diskrétní sinová transformace) a využijeme dopřednou i zpětnou FFT pro odstranění šumu z obrázku (i když se bude jednat o dosti naivní způsob).

První funkcí, se kterou se dnes seznámíme, je funkce nazvaná fft2, kterou lze nalézt v podbalíčku fftpack. Tato funkce je určena pro výpočet FFT dvourozměrného signálu reprezentovaného například dvourozměrnou maticí (NumPy):

fft2(x, shape=None, axes=(-2, -1), overwrite_x=False)
    2-D discrete Fourier transform.
 
    Return the 2-D discrete Fourier transform of the 2-D argument
    `x`.
 
    See Also
    --------
    fftn : for detailed information.
 
    Examples
    --------
    >>> import numpy as np
    >>> from scipy.fftpack import fft2, ifft2
    >>> y = np.mgrid[:5, :5][0]
    >>> y
    array([[0, 0, 0, 0, 0],
           [1, 1, 1, 1, 1],
           [2, 2, 2, 2, 2],
           [3, 3, 3, 3, 3],
           [4, 4, 4, 4, 4]])
    >>> np.allclose(y, ifft2(fft2(y)))

Výpočet 2D FFT pro rastrový obrázek reprezentovaný ve stupních šedi

Jak jsme si již řekli v úvodní kapitole, může být vstupní signál, pro který se 2D FFT počítá, reprezentován maticí. A tato matice může obsahovat hodnoty pixelů rastrového obrázku. Již předminule jsme se seznámili s datovými sadami knihovny SciPy, mj. i s obrázkem schodiště. Připomeňme si, že se jedná o rastrový obrázek o rozměrech 512×512 pixelů, který pochopitelně můžeme podrobit 2D FFT:

import numpy as np
 
import scipy.datasets as datasets
from scipy import fftpack
 
 
# načtení matice
ascent = datasets.ascent()
 
# výpočet 2D FFT
ascent_fft = fftpack.fft2(ascent)
 
np.info(ascent_fft)
 
a = np.abs(ascent_fft)
print(np.min(a))
print(np.max(a))

Po spuštění tohoto skriptu se vypíšou základní informace o výsledku 2D FFT:

class:  ndarray
shape:  (512, 512)
strides:  (8192, 16)
itemsize:  16
aligned:  True
contiguous:  True
fortran:  False
data pointer: 0x7fb586364010
byteorder:  little
byteswap:  False
type: complex128

Z tohoto výpisu je zřejmé, že výsledkem je matice 512×512 komplexních hodnot. Dále se pro zajímavost vypíše maximální a minimální absolutní hodnota získaná z této matice:

4.461849501070654
22932324.0

Zobrazení výsledku výpočtu 2D FFT

Pochopitelně nám nic nebrání v tom, abychom si nechali zobrazit výsledek výpočtu 2D FFT. Víme již, že výsledkem je matice, takže by mělo být možné zobrazit například velikosti komplexních hodnot, popř. fázi atd. Samotná implementace je snadná:

import numpy as np
 
import scipy.datasets as datasets
from scipy import fftpack
 
import matplotlib.pyplot as plt
 
 
# načtení matice
ascent = datasets.ascent()
 
# výpočet 2D FFT
ascent_fft = fftpack.fft2(ascent)
 
# zobrazení výsledku, změna měřítka
plt.imshow(np.abs(ascent_fft))
 
plt.title("2D FFT")
 
# uložení matice do rastrového obrázku
plt.savefig("fft2d_2.png")
 
# zobrazení grafu
plt.show()

Ovšem výsledek nebude příliš oslňující. Proč tomu tak je, si řekneme v navazující kapitole:

Využití logaritmické stupnice při vizualizaci 2D FFT

Obrázek z předchozí kapitoly ve skutečnosti neobsahoval pixely s totožnou barvou; rozdílné hodnoty byly k nalezení ve všech čtyřech rozích (FFT je totiž opět posunuto kvůli záporným frekvencím). Lepšího výsledku dosáhneme tehdy, pokud se použije logaritmické měřítko (škála) při mapování mezi hodnotou 2D FFT a barvou pixelu. Jedno z řešení může být založeno na funkci pro výpočet logaritmu:

import numpy as np
 
import scipy.datasets as datasets
from scipy import fftpack
 
import matplotlib.pyplot as plt
 
 
# načtení matice
ascent = datasets.ascent()
 
# výpočet 2D FFT
ascent_fft = fftpack.fft2(ascent)
 
# zobrazení výsledku, změna měřítka
plt.imshow(np.abs(np.log10(ascent_fft)))
 
plt.title("2D FFT")
 
# uložení matice do rastrového obrázku
plt.savefig("fft2d_3.png")
 
# zobrazení grafu
plt.show()

Nyní bude výsledný obrázek vypadat mnohem lépe:

Vycentrování výsledků FFT (důvod jsme si uvedli již minule):

import numpy as np
 
import scipy.datasets as datasets
from scipy import fftpack
from scipy.fft import fftshift
 
import matplotlib.pyplot as plt
 
 
# načtení matice
ascent = datasets.ascent()
 
# výpočet 2D FFT
ascent_fft = fftpack.fft2(ascent)
 
# posun ve frekvenční oblasti
ascent_fft = fftshift(ascent_fft)
 
# zobrazení výsledku, změna měřítka
plt.imshow(np.abs(np.log10(ascent_fft)))
 
plt.title("2D FFT")
 
# uložení matice do rastrového obrázku
plt.savefig("fft2d_3.png")
 
# zobrazení grafu
plt.show()

Nyní bude výsledek vypadat mnohem použitelněji:

Logaritmické měřítko použité při normalizaci barev

Alternativně je možné použít logaritmické měřítko v průběhu normalizace barev, které je prováděno knihovnou Matplotlib při vizualizaci matice. Povšimněte si dalšího pojmenovaného parametru předaného funkci plt.imshow:

import numpy as np
 
import scipy.datasets as datasets
from scipy import fftpack
 
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
 
 
# načtení matice
ascent = datasets.ascent()
 
# výpočet 2D FFT
ascent_fft = fftpack.fft2(ascent)
 
# zobrazení výsledku
plt.imshow(np.abs(ascent_fft), norm=LogNorm(vmin=5))
plt.colorbar()
 
plt.title("2D FFT")
 
# uložení matice do rastrového obrázku
plt.savefig("fft2d_4.png")
 
# zobrazení grafu
plt.show()

Výsledný obrázek:

Opět bude vhodnější posunout hodnoty ve frekvenční oblasti s využitím funkce fftshift:

import numpy as np
 
import scipy.datasets as datasets
from scipy import fftpack
 
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
 
 
# načtení matice
ascent = datasets.ascent()
 
# výpočet 2D FFT
ascent_fft = fftpack.fft2(ascent)
 
# posun ve frekvenční oblasti
ascent_fft = fftshift(ascent_fft)
 
# zobrazení výsledku
plt.imshow(np.abs(ascent_fft), norm=LogNorm(vmin=5))
plt.colorbar()
 
plt.title("2D FFT")
 
# uložení matice do rastrového obrázku
plt.savefig("fft2d_4B.png")
 
# zobrazení grafu
plt.show()

Nyní bude výsledek vypadat takto:

Zpětná transformace z frekvenční oblasti do oblasti časové

Připomeňme si, že pro jednorozměrné signály je možné zpětnou rychlou Fourierovu transformaci vypočítat s využitím funkce nazvané příznačně ifft. Pro dvourozměrné signály (a tedy i rastrové obrázky) existuje podobná funkce nazvaná ifft2. Tato funkce je opět definována v podbalíčku fftpack:

ifft2(x, shape=None, axes=(-2, -1), overwrite_x=False)
    2-D discrete inverse Fourier transform of real or complex sequence.
 
    Return inverse 2-D discrete Fourier transform of
    arbitrary type sequence x.
 
    See `ifft` for more information.
 
    See Also
    --------
    fft2, ifft
 
    Examples
    --------
    >>> import numpy as np
    >>> from scipy.fftpack import fft2, ifft2
    >>> y = np.mgrid[:5, :5][0]
    >>> y
    array([[0, 0, 0, 0, 0],
           [1, 1, 1, 1, 1],
           [2, 2, 2, 2, 2],
           [3, 3, 3, 3, 3],
           [4, 4, 4, 4, 4]])
    >>> np.allclose(y, fft2(ifft2(y)))

V dalším demonstračním příkladu je nejdříve rastrový obrázek převeden z časové oblasti do oblasti frekvenční funkcí fft2 a posléze je proveden převod zpět, nyní funkcí ifft2. Výsledný obrázek je následně vizualizován, ovšem s tím, že je použita barvová paleta obsahující pouze stupně šedi:

import numpy as np
 
import scipy.datasets as datasets
from scipy import fftpack
 
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
 
 
# načtení matice
ascent = datasets.ascent()
 
# výpočet 2D FFT
ascent_fft = fftpack.fft2(ascent)
 
ascent_time = fftpack.ifft2(ascent_fft)
 
# zobrazení výsledku
plt.imshow(ascent_time.real, cmap="gray")
 
# uložení matice do rastrového obrázku
plt.savefig("fft2d_5.png")
 
# zobrazení grafu
plt.show()

Výsledek by měl (až na malé chyby, které při výpočtu vzniknou) odpovídat originálu:

Využití 2D FFT při filtraci obrazu

Dopřednou a zpětnou rychlou Fourierovu transformaci 2D obrazů je možné využít i pro zpracování obrazů, například pro odstranění šumu, nalezení textu ve fotografiích, nalezení či zvýraznění hran atd. Namísto klasických konvolučních filtrů se v těchto případech obraz nejdříve převede do frekvenční oblasti rychlou FFT, což nám umožní nějakým způsobem manipulovat s jeho spektrálními vlastnostmi. Například můžeme omezit vyšší frekvence (primitivní odstranění šumu), naopak je posílit (detekce hran), nalézt oblasti se specifickými závislostmi mezi frekvencemi v horizontálním a vertikálním směru (detekce písma) atd. Po provedení úprav ve frekvenční oblasti se obraz převede zpět do oblasti časové zpětnou rychlou Fourierovou transformací.

V navazujících kapitolách si ukážeme velmi jednoduchý (a poměrně primitivní) postup pro odstranění šumu z obrazu resp. pro odstranění vyšších frekvencí.

Úprava obrazu ve frekvenční oblasti: odstranění vyšších frekvencí

V předchozí kapitole jsme si naznačili, jakým způsobem je možné provádět filtrace 2D signálu. Ukažme si tedy slíbený jednoduchý filtr, který se pokusí z obrazu odstranit vyšší frekvence. Využijeme toho, že po provedení 2D FFT budou nižší frekvence umístěny okolo rohů (ve frekvenční oblasti). Tyto části nebudeme modifikovat; naopak vynulujeme střední oblasti (ve tvaru kříže, resp. překrývajícího se horizontální a vertikálního pruhu). Na tomto místě je vhodné poznamenat, že se jedná o značně primitivní způsob filtrace, ovšem jde pouze o základní ukázku možnosti modifikace signálu ve frekvenční oblasti:

import numpy as np
 
import scipy.datasets as datasets
from scipy import fftpack
 
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
 
 
# načtení matice
ascent = datasets.ascent()
 
# výpočet 2D FFT
ascent_fft = fftpack.fft2(ascent)
 
high_freq_limit = 100
 
# odstranění vyšších frekvencí
ascent_fft[high_freq_limit:512-high_freq_limit] = 0
ascent_fft[:, high_freq_limit:512-high_freq_limit] = 0
 
ascent_time = fftpack.ifft2(ascent_fft)
 
# zobrazení výsledku
plt.imshow(ascent_time.real, cmap="gray")
 
# uložení matice do rastrového obrázku
plt.savefig("fft2d_6.png")
 
# zobrazení grafu
plt.show()

Rastrový obrázek po výše popsaných úpravách bude vypadat následovně:

Vliv odstranění nižších frekvencí z obrazu

V předchozím skriptu jsme použili konstantu high_freq_limit, jejíž modifikací lze řídit, od jaké frekvence dojde k filtraci. Čím nižší je tato hodnota, tím větší část frekvencí bude vynulována. Pokusme se tedy o vynulování nižších frekvencí:

import numpy as np
 
import scipy.datasets as datasets
from scipy import fftpack
 
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
 
 
# načtení matice
ascent = datasets.ascent()
 
# výpočet 2D FFT
ascent_fft = fftpack.fft2(ascent)
 
high_freq_limit = 50
 
# odstranění vyšších frekvencí
ascent_fft[high_freq_limit:512-high_freq_limit] = 0
ascent_fft[:, high_freq_limit:512-high_freq_limit] = 0
 
ascent_time = fftpack.ifft2(ascent_fft)
 
# zobrazení výsledku
plt.imshow(ascent_time.real, cmap="gray")
 
# uložení matice do rastrového obrázku
plt.savefig("fft2d_7.png")
 
# zobrazení grafu
plt.show()

Ve výsledném obrazu je patrné, že dojde k většímu rozmazání a navíc i ke vzniku artefaktů:

Můžeme se pokusit vynulovat ještě nižší frekvence:

import numpy as np
 
import scipy.datasets as datasets
from scipy import fftpack
 
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
 
 
# načtení matice
ascent = datasets.ascent()
 
# výpočet 2D FFT
ascent_fft = fftpack.fft2(ascent)
 
high_freq_limit = 25
 
# odstranění vyšších frekvencí
ascent_fft[high_freq_limit:512-high_freq_limit] = 0
ascent_fft[:, high_freq_limit:512-high_freq_limit] = 0
 
ascent_time = fftpack.ifft2(ascent_fft)
 
# zobrazení výsledku
plt.imshow(ascent_time.real, cmap="gray")
 
# uložení matice do rastrového obrázku
plt.savefig("fft2d_8.png")
 
# zobrazení grafu
plt.show()

Nyní bude obraz do značné míry znehodnocen, protože se zcela ztratí jakékoli detaily:

Přidání šumu do vstupního obrázku

V dalším kroku do vstupního obrázku vložíme umělý šum. Prozatím použijeme tu nejjednodušší variantu, tj. přidání pseudonáhodných hodnot s předepsaným rozložením ke všem pixelům na obrázku. K tomuto účelu lze použít například funkci numpy.random.normal, které se kromě parametrů náhodných hodnot musí předat i tvar (shape) výsledného n-dimenzionálního pole se šumem. Tvar musí být pochopitelně totožný s tvarem vstupního obrázku, aby bylo možné prvky obou polí následně sečíst:

import numpy as np
 
import scipy.datasets as datasets
import matplotlib.pyplot as plt
 
# načtení matice
ascent = datasets.ascent()
 
# vygenerování šumu
noise = np.random.normal(0, 60, ascent.shape)
 
# zašumění obrázku
ascent = ascent + noise
 
# zobrazení matice
plt.imshow(ascent, cmap="gray")
 
# uložení matice do rastrového obrázku
plt.savefig("fft2d_9.png")
 
# zobrazení grafu
plt.show()

Po výše uvedených úpravách vznikne následující rastrový obrázek:

Pokus o odstranění šumu s využitím FFT s ořezáním vyšších frekvencí

Zašuměný obrázek vytvořený předchozím skriptem převedeme do frekvenční oblasti s využitím 2D FFT, vynulujeme prvky odpovídající vyšším frekvencím a následně provedeme zpětný převod do časové oblasti. Výsledkem by měl být rastrový obrázek s částečně odstraněným šumem (který má vysokou frekvenci):

import numpy as np
 
import scipy.datasets as datasets
from scipy import fftpack
 
import matplotlib.pyplot as plt
 
 
# načtení matice
ascent = datasets.ascent()
 
# vygenerování šumu
noise = np.random.normal(0, 60, ascent.shape)
 
# zašumění obrázku
ascent = ascent + noise
 
# výpočet 2D FFT
ascent_fft = fftpack.fft2(ascent)
 
high_freq_limit = 100
 
ascent_fft[high_freq_limit:512-high_freq_limit] = 0
ascent_fft[:, high_freq_limit:512-high_freq_limit] = 0
 
ascent_time = fftpack.ifft2(ascent_fft)
 
# zobrazení výsledku
plt.imshow(ascent_time.real, cmap="gray")
 
# uložení matice do rastrového obrázku
plt.savefig("fft2d_A.png")
 
# zobrazení grafu
plt.show()

A takto by měl vypadat výsledek (mimochodem nepříliš ucházející):

Odstranění specifických chyb z rastrového obrázku s reálným poškozením

Šum, který byl do obrázku přidán, byl dosti umělý. Zkusme tedy odlišný vstupní obrázek, konkrétně obrázek poškozený mj. i chybami kamery, chybami při přenosu atd. Jedná se o tuto fotografii:

Opět provedeme operace, které již dobře známe: převod do frekvenční oblasti, vynulování vyšších frekvencí a převod zpět:

import numpy as np
 
import scipy.datasets as datasets
from scipy import fftpack
 
import matplotlib.pyplot as plt
 
original = plt.imread("moonlanding.png").astype(float)
 
# výpočet 2D FFT
original_fft = fftpack.fft2(original)
 
high_freq_limit = 100
r, c = original_fft.shape
 
original_fft[high_freq_limit:r-high_freq_limit] = 0
original_fft[:, high_freq_limit:c-high_freq_limit] = 0
 
filtered = fftpack.ifft2(original_fft)
 
# zobrazení výsledku
plt.imshow(filtered.real, cmap="gray")
 
# uložení matice do rastrového obrázku
plt.savefig("fft2d_B.png")
 
# zobrazení grafu
plt.show()

Výsledek vlastně není vůbec špatný, zvláště když si uvědomíme, jak primitivní filtr jsme vlastně použili:

Diskrétní kosinová transformace

Minule jsme se v závěrečné části článku zmínili o diskrétní kosinové transformaci (DCT). Připomeňme si, že tato transformace, která má poměrně úzký vztah k diskrétní Fourierově transformaci (DCT ovšem používá pouze reálné báze), se poměrně často používá při zpracování rastrových obrázků či videa a nalezneme ji taktéž ve formátu JFIF (JPEG). V knihovně SciPy je samozřejmě implementována funkce pro výpočet DCT. Tato funkce se (nepřekvapivě) jmenuje dct a použít ji lze na signál o libovolném množství rozměrů a dokonce je možné si vybrat i typ DCT (I až IV):

dct(x, type=2, n=None, axis=-1, norm=None, overwrite_x=False, workers=None, orthogonalize=None)
    Return the Discrete Cosine Transform of arbitrary type sequence x.
 
    Parameters
    ----------
    x : array_like
        The input array.
    type : {1, 2, 3, 4}, optional
        Type of the DCT (see Notes). Default type is 2.
    n : int, optional
        Length of the transform.  If ``n < x.shape[axis]``, `x` is
        truncated.  If ``n > x.shape[axis]``, `x` is zero-padded. The
        default results in ``n = x.shape[axis]``.
    axis : int, optional
        Axis along which the dct is computed; the default is over the
        last axis (i.e., ``axis=-1``).
    norm : {"backward", "ortho", "forward"}, optional
        Normalization mode (see Notes). Default is "backward".
    overwrite_x : bool, optional
        If True, the contents of `x` can be destroyed; the default is False.
    workers : int, optional
        Maximum number of workers to use for parallel computation. If negative,
        the value wraps around from ``os.cpu_count()``.
        See :func:`~scipy.fft.fft` for more details.
    orthogonalize : bool, optional
        Whether to use the orthogonalized DCT variant (see Notes).
        Defaults to ``True`` when ``norm="ortho"`` and ``False`` otherwise.

Ukázka aplikace DCT na jednoduchý signál reprezentovaný čtyřmi diskrétními hodnotami:

from scipy.fft import dct
 
import numpy as np
 
t = np.array([1, 2, 3, 4])
 
# diskrétní kosinová transformace
f = dct(t)
 
for x in f:
    print(x)

Výsledkem je vektor se čtyřmi reálnými (!) hodnotami:

20.0
-6.308644059797899
0.0
-0.4483415291679651

V knihovně SciPy existuje i funkce nazvaná idct, která počítá zpětnou DCT. Ostatně si můžeme snadno ověřit, že idct(dct(t)) vrátí původní signál t:

from scipy.fft import dct, idct
 
import numpy as np
 
t = np.array([1, 2, 3, 4])
 
f = dct(t)
 
for x in f:
    print(x)
 
print("-" * 10)
 
t2 = idct(f)
 
for x in t2:
    print(x)

Výsledky:

20.0
-6.308644059797899
0.0
-0.4483415291679651
----------
1.0000000000000002
2.0
3.0
4.0

Vizualizace diskrétní kosinové transformace vstupního periodického signálu

Jak vlastně vypadá způsob transformace s využitím DCT, si můžeme ověřit na jednoduchém signálu, který se v časové oblasti skládá ze součtu dvou kosinusovek s různými frekvencemi (frekvence druhé kosinusovky je desetkrát větší, než kosinusovky první). Signál je periodický, ovšem zaznamenáme pouze jednu periodu:

from scipy.fft import dct
 
import numpy as np
import matplotlib.pyplot as plt
 
# signál
t = np.cos(np.linspace(0, 2.0*np.pi, 100)) + 0.5*np.cos(np.linspace(0, 20*np.pi, 100))
 
# vykreslení signálu
plt.plot(t)
 
# uložení grafu s průběhem signálu
plt.savefig("dct_3.png")

Průběh jedné periody signálu:

Nyní si necháme vypočítat DCT tohoto signálu a výsledek zobrazíme (ve frekvenčním spektru):

from scipy.fft import dct
 
import numpy as np
import matplotlib.pyplot as plt
 
# signál
t = np.cos(np.linspace(0, 2.0*np.pi, 100)) + 0.5*np.cos(np.linspace(0, 20*np.pi, 100))
 
# diskrétní kosinová transformace
f = dct(t)
 
# vykreslení výsledného signálu
plt.plot(f)
 
# uložení grafu s průběhem signálu
plt.savefig("dct_4.png")

Výsledek výpočtu DCT:

Diskrétní kosinová transformace signálu se zápornými koeficienty ve frekvenční oblasti

V dalším kroku ve vstupním signálu otočíme druhou kosinusovku (budeme ji odečítat, nikoli přičítat):

from scipy.fft import dct
 
import numpy as np
import matplotlib.pyplot as plt
 
# signál
t = np.cos(np.linspace(0, 2.0*np.pi, 100)) - 0.5*np.cos(np.linspace(0, 20*np.pi, 100))
 
# vykreslení signálu
plt.plot(t)
 
# uložení grafu s průběhem signálu
plt.savefig("dct_5.png")

V časové oblasti si této změny sice můžeme všimnout, ale je dosti nenápadná:

Ovšem v oblasti frekvenční je tomu zcela jinak:

from scipy.fft import dct
 
import numpy as np
import matplotlib.pyplot as plt
 
# signál
t = np.cos(np.linspace(0, 2.0*np.pi, 100)) - 0.5*np.cos(np.linspace(0, 20*np.pi, 100))
 
# diskrétní kosinová transformace
f = dct(t)
 
# vykreslení výsledného signálu
plt.plot(f)
 
# uložení grafu s průběhem signálu
plt.savefig("dct_6.png")

Nyní je jasně patrné, že se signál skládá ze dvou periodických složek, druhá má opačné znaménko (a poloviční velikost):

Diskrétní kosinová transformace 2D obrázku

DCT je možné aplikovat na 2D obrázek, což si pochopitelně ukážeme na demonstračním příkladu. Podobně jako v případě výpočtu FFT, i zde použijeme obrázek schodiště, který je součástí standardních testovacích sad knihovny SciPy. Obrázek běžným způsobem načteme (jedná se o 2D matici), aplikujeme na ní DCT a výsledek přímo zobrazíme:

import numpy as np
 
import scipy.datasets as datasets
from scipy.fft import dct
 
import matplotlib.pyplot as plt
 
 
# načtení matice
ascent = datasets.ascent()
 
# výpočet 2D DCT
ascent_dct = dct(ascent)
 
# zobrazení výsledku, změna měřítka
plt.imshow(ascent_dct)
 
plt.title("2D DCT")
 
# uložení matice do rastrového obrázku
plt.savefig("dct_7.png")
 
# zobrazení grafu
plt.show()

Výsledek ukazuje, že máme podobný problém, jako při výpočtu FFT – lineární mapování mezi hodnotami a barvou:

Musíme tedy použít nelineární mapování, například logaritmickou stupnici. Opět se nejedná o nic nového, protože jsme podobnou operaci prováděli i při vizualizaci výsledků FFT:

import numpy as np
 
import scipy.datasets as datasets
from scipy.fft import dct
 
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
 
 
# načtení matice
ascent = datasets.ascent()
 
# výpočet 2D DCT
ascent_dct = dct(ascent)
 
# zobrazení výsledku, změna měřítka
plt.imshow(ascent_dct, norm=LogNorm(vmin=5))
 
plt.title("2D DCT")
 
# uložení matice do rastrového obrázku
plt.savefig("dct_8.png")
 
# zobrazení grafu
plt.show()

V tomto případě bude výsledek aplikace DCT jasně viditelný:

Jen pro zajímavost (podrobnosti si řekneme příště) se podívejme na DCT obrázku po jeho rozmazání:

import numpy as np
 
import scipy.datasets as datasets
from scipy import ndimage
from scipy.fft import dct
 
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
 
 
# načtení matice
ascent = datasets.ascent()
 
blurred = ndimage.gaussian_filter(ascent, sigma=2.0)
 
# výpočet 2D DCT
blurred_dct = dct(blurred)
 
# zobrazení výsledku, změna měřítka
plt.imshow(blurred_dct, norm=LogNorm(vmin=5))
 
plt.title("2D DCT")
 
# uložení matice do rastrového obrázku
plt.savefig("dct_8B.png")
 
# zobrazení grafu
plt.show()

Diskrétní sinová transformace

V závěrečné části dnešního článku se ještě v krátkosti zmíníme o diskrétní sinové transformaci (DST). Tato transformace má mnoho prvků společných s výše popsanou diskrétní kosinovou transformací, ovšem zatímco DCT používá pouze kosinusové (tedy sudé) báze, pak DST používá pouze báze sinusové (tedy liché). Výsledky DCT lze vhodnými lineárními operacemi převést na DST a naopak. V praxi se setkáme spíše s DCT, která je vhodnější pro operace s obrázky atd. Nicméně se vraťme k DST. Ta je v knihovně SciPy realizována funkcí pojmenovanou dst:

dst(x, type=2, n=None, axis=-1, norm=None, overwrite_x=False, workers=None, orthogonalize=None)
    Return the Discrete Sine Transform of arbitrary type sequence x.
 
    Parameters
    ----------
    x : array_like
        The input array.
    type : {1, 2, 3, 4}, optional
        Type of the DST (see Notes). Default type is 2.
    n : int, optional
        Length of the transform. If ``n < x.shape[axis]``, `x` is
        truncated.  If ``n > x.shape[axis]``, `x` is zero-padded. The
        default results in ``n = x.shape[axis]``.
    axis : int, optional
        Axis along which the dst is computed; the default is over the
        last axis (i.e., ``axis=-1``).
    norm : {"backward", "ortho", "forward"}, optional
        Normalization mode (see Notes). Default is "backward".
    overwrite_x : bool, optional
        If True, the contents of `x` can be destroyed; the default is False.
    workers : int, optional
        Maximum number of workers to use for parallel computation. If negative,
        the value wraps around from ``os.cpu_count()``.
        See :func:`~scipy.fft.fft` for more details.
    orthogonalize : bool, optional
        Whether to use the orthogonalized DST variant (see Notes).
        Defaults to ``True`` when ``norm="ortho"`` and ``False`` otherwise.
 
        .. versionadded:: 1.8.0

Ukažme si výpočet DST na signálu, který jsme již použili ve čtrnácté kapitole při vizualizaci výsledku DCT. Výsledky obou transformací bude možné snadno vizuálně porovnat:

from scipy.fft import dst
 
import numpy as np
import matplotlib.pyplot as plt
 
# signál
t = np.cos(np.linspace(0, 2.0*np.pi, 100)) + 0.5*np.cos(np.linspace(0, 20*np.pi, 100))
 
# diskrétní sinová transformace
f = dst(t)
 
# vykreslení výsledného signálu
plt.plot(f)
 
# uložení grafu s průběhem signálu
plt.savefig("dst_1.png")

Výsledný signál ve frekvenční oblasti bude vypadat následovně:

Diskrétní sinová transformace 2D obrázku

I DST je možné, podobně jako FFT či DCT, aplikovat na dvourozměrné obrázky. Pro otestování této funkcionality opět použijeme testovací datovou sadu s obrázkem schodiště:

import numpy as np
 
import scipy.datasets as datasets
from scipy.fft import dst
 
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
 
 
# načtení matice
ascent = datasets.ascent()
 
# výpočet 2D dst
ascent_dst = dst(ascent)
 
# zobrazení výsledku, změna měřítka
plt.imshow(ascent_dst, norm=LogNorm(vmin=5))
 
plt.title("2D dst")
 
# uložení matice do rastrového obrázku
plt.savefig("dst_3.png")
 
# zobrazení grafu
plt.show()

Podívejme se na vypočtený výsledek (a opět ho můžeme porovnat s výsledkem aplikace DCT):

Repositář s demonstračními příklady

Všechny demonstrační příklady popsané v tomto článku naleznete i v repositáři https://github.com/tisnik/most-popular-python-libs. Následují odkazy na jednotlivé příklady:

Odkazy na Internetu

1. SciPy homepage
   https://scipy.org/
2. SciPy (Wikipedia)
   https://en.wikipedia.org/wiki/SciPy
3. The Hertzsprung–Russell diagram
   https://scipython.com/book2/chapter-9-data-analysis-with-pandas/problems/p92/the-hertzsprung-russell-diagram/
4. Linear algebra (scipy.linalg)
   https://docs.scipy.org/doc/sci­py/reference/linalg.html
5. Frequently Asked Questions – SciPy
   https://scipy.org/faq/
6. SciPy – Introduction
   https://www.tutorialspoin­t.com/scipy/scipy_introduc­tion.htm
7. LAPACK — Linear Algebra PACKage
   https://www.netlib.org/lapack/
8. LAPACK (Wikipedia)
   https://en.wikipedia.org/wiki/LAPACK
9. LAPACK na GitHubu
   https://github.com/Reference-LAPACK/lapack
10. SciPy in Python
    https://pythonguides.com/scipy/
11. scipy.linalg.det
    https://docs.scipy.org/doc/sci­py/reference/generated/sci­py.linalg.det.html#scipy.li­nalg.det
12. scipy.linalg.inv
    https://docs.scipy.org/doc/sci­py/reference/generated/sci­py.linalg.inv.html#scipy.li­nalg.inv
13. scipy.linalg.solve
    https://docs.scipy.org/doc/sci­py/reference/generated/sci­py.linalg.inv.html#scipy.li­nalg.inv
14. Algebra
    https://cs.wikipedia.org/wiki/Algebra
15. Lineární algebra
    https://cs.wikipedia.org/wi­ki/Line%C3%A1rn%C3%AD_alge­bra
16. Lineární rovnice
    https://cs.wikipedia.org/wi­ki/Line%C3%A1rn%C3%AD_rov­nice
17. Soustava lineárních rovnic
    https://cs.wikipedia.org/wi­ki/Soustava_line%C3%A1rn%C3%AD­ch_rovnic
18. Norma matice
    https://cs.wikipedia.org/wi­ki/Norma_matice
19. Matrix norm
    https://en.wikipedia.org/wi­ki/Matrix_norm
20. Norma (vektoru)
    https://cs.wikipedia.org/wi­ki/Norma_(matematika)
21. Frobeniův skalární součin
    https://cs.wikipedia.org/wi­ki/Frobeni%C5%AFv_skal%C3%A1rn%C3%AD_sou­%C4%8Din
22. BLAS (Basic Linear Algebra Subprograms)
    https://www.netlib.org/blas/
23. Basic Linear Algebra Subprograms (Wikipedia)
    https://en.wikipedia.org/wi­ki/Basic_Linear_Algebra_Sub­programs
24. Operace s daty uloženými v binárních souborech v knihovnách NumPy a Pandas
    https://www.root.cz/clanky/operace-s-daty-ulozenymi-v-binarnich-souborech-v-knihovnach-numpy-a-pandas/
25. Operace s daty uloženými v binárních souborech v knihovnách NumPy a Pandas (dokončení)
    https://www.root.cz/clanky/operace-s-daty-ulozenymi-v-binarnich-souborech-v-knihovnach-numpy-a-pandas-dokonceni/
26. NumPy Home Page
    http://www.numpy.org/
27. NumPy v1.10 Manual
    http://docs.scipy.org/doc/num­py/index.html
28. NumPy (Wikipedia)
    https://en.wikipedia.org/wiki/NumPy
29. OpenBLAS: An optimized BLAS library
    https://www.openblas.net/
30. Integrovaná vývojová prostředí ve Fedoře: IPython a IPython Notebook
    http://mojefedora.cz/integrovana-vyvojova-prostredi-ve-fedore-ipython-a-ipython-notebook/
31. Integrovaná vývojová prostředí ve Fedoře: praktické použití IPython Notebooku a knihovny Numpy
    http://mojefedora.cz/integrovana-vyvojova-prostredi-ve-fedore-prakticke-pouziti-ipython-notebooku-a-knihovny-numpy/
32. Integrovaná vývojová prostředí ve Fedoře: praktické použití IPython Notebooku a knihovny Numpy (2.část)
    http://mojefedora.cz/integrovana-vyvojova-prostredi-ve-fedore-prakticke-pouziti-ipython-notebooku-a-knihovny-numpy-2-cast/
33. Integrovaná vývojová prostředí ve Fedoře: vykreslování grafů s využitím knihoven Numpy a matplotlib
    http://mojefedora.cz/integrovana-vyvojova-prostredi-ve-fedore-vykreslovani-grafu-s-vyuzitim-knihoven-numpy-a-matplotlib/
34. Integrovaná vývojová prostředí ve Fedoře: vykreslování grafů s využitím knihoven Numpy a matplotlib (2.část)
    http://mojefedora.cz/integrovana-vyvojova-prostredi-ve-fedore-vykreslovani-grafu-s-vyuzitim-knihoven-numpy-a-matplotlib-2-cast/
35. Piecewise linear interpolation
    https://docs.scipy.org/doc//sci­py/tutorial/interpolate/1D­.html
36. Statistical functions (scipy.stats)
    https://docs.scipy.org/doc/sci­py/reference/stats.html
37. scipy.stats.linregress
    https://docs.scipy.org/doc/sci­py/reference/generated/sci­py.stats.linregress.html
38. numpy.poly1d
    https://numpy.org/doc/sta­ble/reference/generated/num­py.poly1d.html
39. scipy.optimize.curve_fit
    https://docs.scipy.org/doc/sci­py/reference/generated/sci­py.optimize.curve_fit.html#sci­py.optimize.curve_fit
40. scipy.interpolate.make_smoothing_spline
    https://docs.scipy.org/doc/sci­py/reference/generated/sci­py.interpolate.make_smoot­hing_spline.html#scipy.in­terpolate.make_smoothing_spli­ne
41. Famous Curves Index
    https://mathshistory.st-andrews.ac.uk/Curves/
42. Curve (Wikipedia)
    https://en.wikipedia.org/wiki/Curve
43. Mathematical curves
    https://www.2dcurves.com/index.html
44. Curves (Wolfram MathWorld)
    https://mathworld.wolfram­.com/topics/Curves.html
45. Smooth Curve (Wolfram MathWorld)
    https://mathworld.wolfram­.com/SmoothCurve.html
46. Spirals (Wolfram MathWorld)
    https://mathworld.wolfram­.com/topics/Spirals.html
47. An Interactive Introduction to Splines
    https://ibiblio.org/e-notes/Splines/Intro.htm
48. Algebraic curve
    https://en.wikipedia.org/wi­ki/Algebraic_curve
49. Transcendental curve
    https://en.wikipedia.org/wi­ki/Transcendental_curve
50. Algebraic Curves
    https://mathworld.wolfram­.com/topics/AlgebraicCurves­.html
51. Elliptic Curves
    https://mathworld.wolfram­.com/topics/EllipticCurves­.html
52. Curves we (mostly) don't learn in high school (and applications)
    https://www.youtube.com/wat­ch?v=3izFMB91K_Q
53. Catenary arch
    https://en.wikipedia.org/wi­ki/Catenary_arch
54. Parabolic arch
    https://en.wikipedia.org/wi­ki/Parabolic_arch
55. Wattova křivka
    https://www.geogebra.org/m/gNh4bW9r
56. Fifty Famous Curves, Lots of Calculus Questions, And a Few Answers
    https://elepa.files.wordpres­s.com/2013/11/fifty-famous-curves.pdf
57. Faux, I.D. a Pratt, M.J.: Computational Geometry for Design and Manufacture,
    Ellis Horwood Ltd., Wiley & Sons, 1979
58. Wallace A.: Differential Topology,
    Benjamin/Cummings Co., Reading, Massachussetts, USA, 1968
59. Glossary of Bridge Terminology
    http://sdrc.lib.uiowa.edu/en­g/bridges/WaddellGlossary/Glos­sC.htm
60. Brachistochrona
    https://cs.wikipedia.org/wi­ki/Brachistochrona
61. Missions: Cassini
    https://solarsystem.nasa.gov/mis­sions/cassini/overview/
62. Giovanni Domenico Cassini
    https://en.wikipedia.org/wi­ki/Giovanni_Domenico_Cassi­ni
63. Cassini Ovals
    https://mathworld.wolfram­.com/CassiniOvals.html
64. Geocentrismus
    https://cs.wikipedia.org/wi­ki/Geocentrismus
65. Who was Giovanni Cassini?
    https://www.universetoday­.com/130823/who-was-giovanni-cassini/
66. Special plane curves
    http://xahlee.info/Special­PlaneCurves_dir/ConicSecti­ons_dir/conicSections.html
67. Interpolace
    https://mathonline.fme.vut­br.cz/pg/Algoritmy/05_APROX_KRIV­KY.htm
68. Lagrange Polynomial Interpolation
    https://pythonnumericalmet­hods.berkeley.edu/notebook­s/chapter17.04-Lagrange-Polynomial-Interpolation.html
69. Python Program for Lagrange Interpolation Method (with Output)
    https://www.codesansar.com/numerical-methods/python-program-lagrange-interpolation-method.htm
70. Smooth Paths Using Catmull-Rom Splines
    https://qroph.github.io/2018/07/30/smo­oth-paths-using-catmull-rom-splines.html
71. Lecture 11: Linear Interpolation Again – Bézier Curves
    http://www.math.kent.edu/~re­ichel/courses/intr.num.com­p.1/fall09/lecture12/bez.pdf
72. Geometrie/Úvod do křivek
    https://cs.wikibooks.org/wi­ki/Geometrie/%C3%9Avod_do_k%C5%99i­vek
73. B-Spline Curves and Surfaces (1)
    http://www.cad.zju.edu.cn/ho­me/zhx/GM/006/00-bscs1.pdf
74. Praktické ukázky možností aplikace Mandelbulber při tvorbě animací
    https://www.root.cz/clanky/prakticke-ukazky-moznosti-aplikace-mandelbulber-pri-tvorbe-animaci/
75. Kochanek–Bartels spline
    https://en.wikipedia.org/wi­ki/Kochanek%E2%80%93Bartel­s_spline
76. class KochanekBartels
    https://splines.readthedoc­s.io/en/latest/_modules/spli­nes.html#KochanekBartels
77. Konvoluce
    https://cs.wikipedia.org/wi­ki/Konvoluce
78. Korelace
    https://cs.wikipedia.org/wi­ki/Korelace
79. Fourierova transformace
    https://cs.wikipedia.org/wi­ki/Fourierova_transformace
80. Lanczos resampling
    https://en.wikipedia.org/wi­ki/Lanczos_resampling
81. Obrázek se schodištěm
    https://pixnio.com/people/accent-to-the-top
82. Obrázek s mývalem
    https://pixnio.com/fauna-animals/raccoons/raccoon-procyon-lotor
83. Sample-rate conversion
    https://en.wikipedia.org/wiki/Sample-rate_conversion
84. Seriál Grafické formáty
    https://www.root.cz/serialy/graficke-formaty/
85. JPEG – král rastrových grafických formátů?
    https://www.root.cz/clanky/jpeg-kral-rastrovych-grafickych-formatu/
86. Ztrátová komprese obrazových dat pomocí JPEG
    https://www.root.cz/clanky/ztratova-komprese-obrazovych-dat-pomoci-jpeg/
87. Programujeme JPEG: transformace a podvzorkování barev
    https://www.root.cz/clanky/pro­gramujeme-jpeg-transformace-a-podvzorkovani-barev/
88. Programujeme JPEG: diskrétní kosinová transformace (DCT)
    https://www.root.cz/clanky/pro­gramujeme-jpeg-diskretni-kosinova-transformace-dct/
89. Programujeme JPEG: Kvantizace DCT koeficientů
    https://www.root.cz/clanky/pro­gramujeme-jpeg-kvantizace-dct-koeficientu/
90. Programujeme JPEG: Huffmanovo kódování kvantovaných DCT složek
    https://www.root.cz/clanky/pro­gramujeme-jpeg-huffmanovo-kodovani-kvantovanych-dct-slozek/
91. Programujeme JPEG: Interní struktura souborů typu JFIF/JPEG
    https://www.root.cz/clanky/pro­gramujeme-jpeg-interni-struktura-souboru-typu-jfifjpeg/
92. Programujeme JPEG: Načtení informací ze souborů typu JFIF/JPEG
    https://www.root.cz/clanky/pro­gramujeme-jpeg-nacteni-informaci-ze-souboru-typu-jfifjpeg/
93. Programujeme JPEG: Progresivní JPEG a informace EXIF
    https://www.root.cz/clanky/pro­gramujeme-jpeg-progresivni-jpeg-a-informace-exif/
94. The Fourier Analysis –The Fast Fourier Transform (FFT) Method
    https://www.electronics-lab.com/article/the-fourier-analysis-the-fast-fourier-transform-fft-method/
95. Understanding the output of FFT
    https://howthefouriertran­sformworks.com/understanding-the-output-of-an-fft/
96. What factors go into choosing DCT type 1 over type 2 or type 3 etc?
    https://dsp.stackexchange­.com/questions/31611/what-factors-go-into-choosing-dct-type-1-over-type-2-or-type-3-etc
97. The Discrete Cosine Transform (DCT)
    https://ccrma.stanford.edu/~jos/mdft/Dis­crete_Cosine_Transform_DCT­.html
98. Eulerův vzorec
    https://cs.wikipedia.org/wi­ki/Euler%C5%AFv_vzorec
99. Komplexní analýza
    https://cs.wikipedia.org/wi­ki/Komplexn%C3%AD_anal%C3%BDza
100. What is the difference between a Fourier transform and a cosine transform?
     https://dsp.stackexchange­.com/questions/13/what-is-the-difference-between-a-fourier-transform-and-a-cosine-transform
101. Two-dimensional Discrete Cosine Transform as a Linear Transformation
     https://fairyonice.github.io/2D-DCT.html
102. fft2 (Matlab)
     https://www.mathworks.com/hel­p/matlab/ref/fft2.html
103. Image denoising by FFT
     https://scipy-lectures.org/intro/scipy/au­to_examples/solutions/plot_fft_i­mage_denoise.html
104. Discrete sine transform
     https://en.wikipedia.org/wi­ki/Discrete_sine_transform
105. Diskrétní sinová transformace: demo
     https://www.desmos.com/cal­culator/k5jlr0ykzw