FIT ČVUT Course Pages
BI-JUL — Programování v jazyku Julia
Jdi na navigaci předmětu

07: Vizualizace a manipulace s grafikou

Tento notebook je výukovým materiálem v předmětu BI-JUL.21 vyučovaném v zimním semestru akademického roku 2023/2024 Tomášem Kalvodou. Tvorba těchto materiálů byla podpořena NVS FIT.

Hlavní stránkou předmětu, kde jsou i další notebooky a zajímavé informace, je jeho Course Pages stránka.

versioninfo()

1. Vizualizace

Nejen v Julia existuje celá řada nástrojů/knihoven pro vytváření různých grafů, diagramů a vizualizací. V Python ekosystému je asi nejpoúžívanější Matplotlib, který je k dispozici i v Julia (viz PyPlot.jl níže).

Všechny tyto nástroje se snaží dosáhnout podobných cílů a nabízejí různě kvalitní implementace. Dost často bývá otázkou vkusu, který z nástrojů použít. V tomto notebooku ukážeme tři možnosti:

Paleta možností je ale širší. Zmiňme alespoň následující (pro zájemce k prozkoumání):

  • Plots: nadstavba nad několika grafickými balíčky.
  • Gadfly: atraktivní a interaktivní výstup.
  • Winston: malý a rychlý, terminologií vychází s MATLABu.

Možností je opravdu mnoho a tento notebook je spíše experimentální.

Při používání tohoto notebooku doporučuji vždy mezi sekcemi (různými balíčky) restartovat jádro. Moduly mají podobné zaměření a mohlo by dojít ke kolizím pojmenování metod.

1.1 GR.jl

Balíček GR.jl poskytuje Julia rozhranní k GR frameworku. GR byl vyvinut skupinou Scientific IT-Systems z Peter Grünberg Institutu na Forschungszentrum Jülich.

  • Mezi jeho přednosti patří rychlost.
  • Mezi nevýhody bych zařadil nekvalitní dokumentaci a nekonzistenci API. (Chyba může být ale na straně Julia balíčku...)

Tento balíček pravděpodobně nemáte nainstalovaný, takže nejprve tento musíme nainstalovat (případně viz pokyny k instalaci systémových závislostí):

(@v1.6) pkg> add GR

Zde v notebooku projdeme jen ty nejužitečnější části balíčku. Také tento balíček použijeme jako jakéhosi průvodci po možnostech, u ostatních se totiž dost opakují. Zvídavého čtenáře odkazujeme na dokumentaci s ukázkami dalších typů grafů.

using GR

Tvorbu jednotlivých grafů můžeme oddělit zavoláním figure, případně subplot. K dispozici také je savefig.

plot a oplot

Pomocí plot nepřekvapivě vykreslíme 2D graf funkce. Poziční argumenty mohou mít následující význam:

  • stejně dlouhé vektory x a y reprezentující souřadnice bodů,
  • vektor x a funkce pro výpočet yy-nové hodnoty,
  • pouze vektor y, nezávisle proměnná pak odpovídá indexům.

Dále plot přijímá několik keyword argumentů, nejzajímavější asi jsou:

  • title: titulek grafu (alternativně metoda title),
  • xlabel: titulek grafu (alternativně metoda xlabel),
  • ylabel: titulek grafu (alternativně metoda ylabel),
  • xlim, ylim: rozsah jednotlivých os (alternativně metody xlim a ylim),

a metoda

  • legend: popisky jednotlivých grafů.

Následuje několik jednoduchých ukázek:

xs = LinRange(-10, 10, 200)
ys = sin.(xs)

figure()
plot(xs, ys, title="\$y=\\sin(x)\$", xlabel="\$x\$", ylabel="\$y\$", ylim=(-1.2, 1.2))

Pokud chceme vykreslit více dat do jednoho obrázku, použijeme k tomu oplot, nastavení popisků atp. se bohužel bere z posledního volání oplus (lze to obejít nastavením těchto parametrů pomocí metod, ne pomocí keyword hodnot).

legend("\$\\sin\$", "\$\\cos\$")

oplot(xs, cos, title="Trigonometrické funkce", xlabel="\$x\$", ylabel="\$y\$", ylim=(-1.2, 1.2))

Styl čáry by mělo jít kontrolovat pomocí klíčů:

  • setlinetype: druh čáry (viz GR.LINETYPE_...),
  • setlinewidth: šířka čáry,
  • setlinecolorind: barva čáry (???).
figure()

setlinetype(GR.LINETYPE_DOUBLE_DOT)
plot(LinRange(-1, 1, 50), exp, linewidth=4)

Zde jsem se vzdal a proto:

Cvičení: ???

  1. Pokuste* se změnit barvu čáry například na červenou.

* tj. úspěch není nutný.

figure()

plot(LinRange(-1, 1, 50), exp, "--r", linewidth=4)

scatter

Scatter plot vykresluje pouze jednotlivé body a neinterpoluje (lineárně) mezi nimi:

figure()

title("\$a_n = (-1)^n\$")
legend("\$a_n\$")
xlabel("\$n\$")
ylabel("\$a_n\$")

plot([(-1)^n for n=0:9])
figure()

title("\$a_n = (-1)^n\$")
legend("\$a_n\$")
xlabel("\$n\$")
ylabel("\$a_n\$")

x = Float64[j for j=0:9]
y = Float64[(-1)^n for n=0:9]
s = 2*ones(10)
c = [1, 1, 1, 1, 1, 1, 1, 1, 1, 2]

scatter(x, y, s, c, xlim=(-0.2, 9.2), ylim=(-1.2, 1.2))
x
y
s
c
figure()
x = [z/10 for z=0:10]
y = LinRange(0, 1, 11)
s = LinRange(1, 3, 11)
c = LinRange(0, 255, 11)
scatter(x, y, s, c)
x
y
s
c

plot3

plot3 je 3D analog čárového grafu. Přijímá tedy tři vektory x, y a z, jejichž složky popořadě tvoří jednotlivé body křivky.

ts = LinRange(0, 20, 400)
xs = [ sin(t) for t in ts ]
ys = [ cos(t) for t in ts ]
zs = ts

figure()
plot3(xs, ys, zs, xlabel="x", ylabel="y", zlabel="z", title="Šroubovice")

histogram a hexbin, Cvičení

Histogram bere vektor x hodnot a jako keyword argument nbins udávající počet přihrádek. Budeme ho ilustrovat na náhodném generátoru bodů na kružnici (vzpomeňte dřívější diskuzi).

xs = rand(5000) .- 0.5
ys = rand(5000) .- 0.5;
figure()
scatter(xs, ys, xlim=(-1.2,1.2), ylim=(-1.2,1.2))
for j=1:length(xs)
    r = sqrt(xs[j]^2 + ys[j]^2)
    xs[j] /= r
    ys[j] /= r
end
figure(width=400, height=400) # aspect ratio...?
scatter(xs, ys)
φs = [ atan(ys[j], xs[j]) for j=1:length(xs) ];
figure()
histogram(φs, nbins=50, xlabel="\$\\varphi\$", ylabel="hits")
savefig("test.pdf")
savefig("test.png")

Takovýto generátor evidentně preferuje směry odpovídající osám kvadrantů.

Dvourozměrný ekvivalent histogramu (tj. pro dvousložková data) je metoda hexbin, která funguje analogicky, jen rovinu rozděluje na šestiúhelníky.

xs = randn(10^6)
ys = randn(10^6)

figure()
hexbin(xs, ys)

contour a contourf

Contour plot je vizualizace funkce dvou proměnných, která v jejím definičním oboru znázorňuje křivky, kde je daná funkce konstantní (tzv. kontury). Varianta s f v názvu vybarvuje plochy mezi konturami (fill).

Metody očekávají buď trojici vektorů x, y, z, kde x a y definují souřadnice mřížky a z obsahuje funkční hodnoty na této mřížce (tj. z je matice). Místo z můžeme předat funkci dvou proměnných.

figure()
contour(LinRange(-2, 2, 50), LinRange(-1, 1, 50), (x, y) -> 4x^2 + y^2)

Pozor, jak se GR chová k pořadí os!

figure()

xs = LinRange(-10, 10, 100)
ys = LinRange(-15, 15, 100)
zs = [ sin(x)*cos(y) for x in xs, y in ys ]

contourf(xs, ys, zs)
zs

heatmap

heatmap funguje v podstatě stejně jako contour, jen se nesnaží nalézt křivky konstantní hodnoty funkce (což je netriviální! Odpovídající algoritmus marching cubes algorithm (zde squares -- 2D varianta) byl ještě do relativně nedávná doby pod patentem.)

figure()

xs = LinRange(-10, 10, 200)
ys = LinRange(-10, 10, 200)
zs = [ sin(x)*cos(y/3) for x in xs, y in ys ]

heatmap(xs, ys, zs, colormap=GR.COLORMAP_FLAME)
figure()

# Create example data
x = LinRange(-2, 2, 40)
y = LinRange(0, pi, 20)
z = sin.(x') .+ cos.(y)
# Draw the heatmap
heatmap(z)

imshow

Tato metoda je velmi užitečná. Umožňuje nám graficky znázornit matici.

h = [ i + j - 1 for i=1:10, j=1:10 ]

figure()
imshow(h, colormap=GR.COLORMAP_INFERNO)
figure()

xs = LinRange(-10, 10, 200)
ys = LinRange(-10, 10, 200)
zs = [ sin(x)*cos(y/3) for x in xs, y in ys ]

imshow(zs, colormap=GR.COLORMAP_WINTER)

Cvičení: Julia v Julia

Uvažme zobrazení f:C→Cf: \mathbb{C} \to \mathbb{C} zadané předpisem fc(z)=z2+cf_c(z) = z^2 + c, pro pevně zadanou konstantu c∈Cc \in \mathbb{C}. Zkoumáme pak množinu bodů z∈Cz\in\mathbb{C}, pro které je posloupnost (fcn(z))n=1∞(f^n_c(z))_{n=1}^\infty omezená (mocnina označuje opakované skládání zobrazení).

Přibližně se lze k této množině blížit takto: zvolme dostatečně velké R>0R > 0 a maximální počet iterací M∈NM \in \mathbb{N} a pro zz z nějaké zvolené množiny (typicky mřížka) naměřme, kolik iterací musíme provést, než se dostaneme dále než RR od počátku.

Proveďte tento výpočet pro c=−0.8+0.156ic = -0.8 + 0.156i a c=(φ−2)+(φ−1)ic = (\varphi - 2) + (\varphi - 1)i.

function julia(c::Complex{Float64}, R::Float64, reals, imags; imax::Int64=1000)
    data = fill(imax, length(reals), length(imags))
    f(u) = u^2 + c   
    
    for j = axes(reals, 1), k = axes(imags, 1)
        z = reals[j] + 1im * imags[k]
        
        for n = 1:imax
            z = f(z)
            if abs(z) > R
                data[j, k] = n
                break
            end
        end
    end
    
    return transpose(data)
end
res = LinRange(-2, 2, 20)
ims = LinRange(-1, 1, 10)
data1 = julia(-0.8 + 0.156im, 4.0, res, ims)
imshow(data1, colormap=GR.COLORMAP_BLUESCALE)
imshow([0 1 0; 0 0 1])
res = LinRange(-2, 2, 600)
ims = LinRange(-1, 1, 400)
data2 = julia(-0.8 + 0.156im, 2.0, res, ims);
imshow(data2, colormap=GR.COLORMAP_BLUESCALE)
res = LinRange(-0.3, -0.1, 800)
ims = LinRange(-0.35, -0.15, 800)
data2 = julia(-0.8 + 0.156im, 4.0, res, ims);
imshow(data2, colormap=GR.COLORMAP_INFERNO)
using PyPlot
res = LinRange(-2, 2, 600)
ims = LinRange(-1, 1, 400)
data2 = julia(-0.8 + 0.156im, 4.0, res, ims);

fig, ax = plt.subplots()
ax.imshow(data2, extent=[res[1], res[end], ims[1], ims[end]])
res = LinRange(0, 2, 600)
ims = LinRange(0, 1, 400)
data2 = julia(-0.8 + 0.156im, 4.0, res, ims);

fig, ax = plt.subplots()
ax.imshow(data2, extent=[res[1], res[end], ims[1], ims[end]])
res = LinRange(-0.6, 0, 1000)
ims = LinRange(0, 0.6, 1000)
c = (MathConstants.golden - 2) + (MathConstants.golden - 1)*1im
data3 = julia(c, 5.0, res, ims);
fig, ax = plt.subplots()
ax.imshow(data3, extent=[res[1], res[end], ims[1], ims[end]], cmap="Reds")

1.2 PyPlot.jl

PyPlot.jl je v podstatě Julia rozhraní k Matplotlib Pyhon knihovně, přesněji k funkcím v matplotlib.pyplot Python modulu.

Pokud znáte tuto knihovnu, pak by pro vás tato varianta měla být asi nejpříjemnější na používání. Pro ostatní mohou být užitečné tyto taháky.

  • Výhody: lepší dokumentace, obrovská základna uživatelů, ozkoušený a robustní software.
  • Nevýhody: pomalejší, závislost na Pythonu.

Budete potřebovat mít nainstalovánu tuto knihovnu ($ pip install matplotlib) a poté nainstalovat Julia balíček PyPlot.jl:

(@v1.6) pkg> add PyPlot

Pro porovnání ukažme jen základní typ grafu a pak zkusme vytvořit komplikovanější diagram v cvičení.

# Doporučuji restartovat jádro!
using PyPlot

plot

Logika této metody je velmi podobná jako analogické metody v GR.jl výše.

xs = LinRange(0, 10, 200)
ys = [ sin(x) - sin(2x)/2 + sin(3x)/3 - sin(4x)/4 for x in xs ]

f = figure()
plot(xs, ys, color="red", linewidth=4.0, linestyle="--")
grid()
title("Můj úžasný graf")
xlabel("\$x\$")
ylabel("\$y\$");

Výsledné obrázky lze opět snadno uložit.

f.savefig("matplotlib_fig.png")

K dispozici také je několik přednastavených stylů, například:

f = figure()
plt.style.use("seaborn")
plot(xs, ys, color="red", linewidth=4.0, linestyle="--")
title("Můj úžasný graf")
xlabel(L"$x$")
ylabel(L"$y$");

Cvičení: Řešení ODE

V předchozí lekci jsme si ukázali, jak vyřešit obyčejnou diferenciální rovnici y′=Ayy' = \mathbb{A} y, y(0)=y0∈Cny(0) = y_0 \in \mathbb{C}^n, pomocí maticové exponenciály. Řešením byla vektorová funkce y(t)=exp⁡(tA)y0y(t) = \exp(t\mathbb{A}) y_0.

Konkrétně jsme zvolili

A=(0−110).\mathbb{A} = \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix}.

Pokuste se vizualizovat tato řešení (v rovině) pro pár počátečních podmínek.

Poznámka: Viz předchozí notebook.

A = [0 -1; 1 0]

function solution(matrix, initial_condition, timestamps)
    # ...
end
ts = LinRange(0, 3, 100)
s1 = solution(A, [1, 0], ts)
s2 = solution(A, [-0.5, 0], ts)
s3 = solution(A, [1, 1], ts);
f = figure()

scatter(s1[:,1], s1[:,2], color="red")
scatter(s2[:,1], s2[:,2], color="green")
scatter(s3[:,1], s3[:,2], color="blue")
grid()
title("Trajektorie")
xlabel("\$x\$")
ylabel("\$y\$");

1.3 PGFPlotsX.jl

Tento balíček využívá pro generování grafů LaTeX balíček PGFPlots.

  • Výhody: zřejmě nejkvalitnější výstup, vhodný i do publikací.
  • Nevýhody: ideálně vyžaduje familiaritu s PGFPlots, náročnější instalace.

Musíte mít lokálně nainstalovaný LaTeX s balíčkem PGFPlots. Poté jen přidáme Julia balíček

(@v1.9) pkg> add PGFPlotsX
using PGFPlotsX

Takřka nutností pro používání tohoto balíčku je prostudování si PGFPlots manuálu, který je velmi podrobný a názorný.

xs = LinRange(0, 20, 100)
ys = [ cos(x)/sqrt(x+1) for x in xs ];

Makro @pgf nám umožňuje parsovat výrazy, které by nebyly validní Julia kód, například parametry grafu:

options = @pgf {
    title => "Tlumená oscilace",
    xlabel => "\$x\$",
    ylabel => "\$y\$",
    grid => "major"
}

Grafické objekty pak vytvoříme pomocí konstruktorů, jejichž jména kopírují názvy LaTeX maker a prostředí.

g = Axis(options, Plot(Table(xs, ys)))

Můžeme si i prohlédnout LaTeX kód, který se vygeneroval:

print_tex(g)

PGFPlots nabízí bohaté možnosti změny stylu a prezentace:

@pgf Axis(options, Plot({"red", "line width" => "2pt" }, Table(xs, ys)))

PGFPlots zvládá i 3D grafy (v LaTeXu!!!).

@pgf Axis({"colorbar"},
    Plot3({"surf", "samples" => "15", "domain" => "0:1", "y domain" => "-1:1"}, Expression("x^2 - y^2"))
)

1.4 Makie.jl

Do čtvrtice ještě jedna poměrně zajímavá možnost. Pro podrobnější prozkoumání doporučuji vaší pozornosti galerii.

  • Výhody: poměrně velmi dobrá dokumentace se spoustou příkladů.
  • Nevýhody: vleče se při prvním startu.

Balíček nabízí několik backendů: CairoMakie.jl (2D), GLMakie.jl a WGLMakie.jl. Instalace je opět standardní a nebudu ji zde opakovat.

using CairoMakie
xs = LinRange(0, 20, 200)
ys = sin.(xs) .* xs;
fig, ax, plt = lines(xs, ys)
lines!(xs, cos.(xs))

ax.xlabel = "x"
ax.ylabel = "y"
ax.title = "Pokusný graf"

current_figure()

3D grafy

Makie poskytuje dva backendy schopné vykreslovat 3D grafiku, GLMakie.jl a WGLMakie.jl. 

using GLMakie

xs = LinRange(0, 10, 100)
ys = LinRange(0, 15, 100)
zs = [cos(x) * sin(2*y) + x + y for x in xs, y in ys]

surface(xs, ys, zs, axis=(type=Axis3,))

Cvičení: Vizualizace PRNG založeného na logistickém zobrazení

function my_prng(n, seed)
    data = zeros(n)
    data[1] = seed
    
    for j = 2:n
        data[j] = 4data[j-1]*(1 - data[j-1])
    end
    
    return data
end
data = my_prng(10^6, 0.7500000001);

Jak tato posloupnost vypadá?

using CairoMakie
fig, ax, plt = scatter(1:1000, data[1:1000])

ax.xlabel = "n"

current_figure()
hist(data, bins=50, normalization = :pdf)
f(x) = acos(1-2x)/pi
udata = f.(data);
hist(udata, bins=50, normalization = :pdf)

1.5 Interact.jl

Tento balíček umožňuje vytvářet interaktivní vizualizace. Pokyny k instalaci jsou tentokrát malinko komplikovanější, viz Getting Started.

using Interact, CairoMakie
xs = LinRange(-5, 20, 200);
@manipulate for a=0:0.01:1, ω=0:0.01:2
    ys  = [ a * sin(ω * x) for x in xs  ]
    vbox(lines(xs, ys))
end

2. Images.jl, Colors.jl, ImageMagick.jl, ImageMetadata.jl...

Tyto balíčky se specializují přímo na práci s obrázky, jejich vytváření a jejich editaci. Nespadají tedy přímo do kategorie předchozích balíčků. Využijete je třeba při rychlém generování obrázků z matic nebo při strojovém před/zpracování obrázků (viz ImageFeatures.jl, ImageTransformation.jl nebo ImageSegmentation.jl, ImageDraw.jl).

Instalace je opět standardní (]add Images, ]add Colors).

using Images, Colors

2.1 Vytvoření obrázku

img_data = rand(400, 600)
img_data[1:50, 1:50] .= 1
img_data
img = Gray.(img_data)
display(MIME("text/plain"), img)
dump(img)

Obrázky můžeme i snadno ukládat v různých formátech.

using FileIO
save("test.jpg", img)

Nahrávání obrázků

using FileIO, ImageMagick, ImageIO
img_Gauss = load("../homeworks/B211/files/gauss_portret.jpg")
typeof(img_Gauss)
dump(img_Gauss)

Všimněte si speciálního datového typu:

?N0f8

Jeden pixel je barevně reprezentován trojicí čísel (RGB).

p = img_Gauss[1, 1]
typeof(p)
p.r
p.g
p.b
img_Gauss[10, 10]

Například můžeme extrahovat červený kanál a převést ho do datového typu, v kterém se lépe počítá:

imgR = map(x -> Float64(x.r), img_Gauss)
Gray.(imgR)
Gray.([imgR[1:180, :]; imgR[225:end, :]])

3. Cairo.jl a Luxor.jl

Cairo.jl představuje rozhranní k 2D grafické knihovně Cairo. S pomocí této knihovny jsme schopni deklarativním způsobem vytvářet zcela vlastní obrázky, ovšem operuje na poměrně nízké úrovni. Uživatelsky příjemnější je balíček Luxor.jl. Tento balíček nám umožňuje vytvářet rasterovou (PNG) i vektorovou (SVG) grafiku. Jeho výhodou je i poměrně podrobná dokumentace s tutoriály.

using Luxor
@png begin
    text("Toto je můj první obrázek!", Point(0, 0), angle = pi/4, halign = :center)
    rect(-50, -50, 100, 100, action = :stroke)
    circle(Point(0, 0), 50, action = :stroke)
end
@svg begin
    circle(Point(-50, 0), 50, action=:fill)
    circle(Point(50, 0), 60, action=:fill)
    rulers()
end

Pozor, souřadný systém je orientován v grafice obvyklým způsobem "hlavou dolů".

Vedle výše použitých primitiv rect a circ máme k dispozici spoustu dalších nástrojů.

@svg begin
    background("blue")
    sethue("red")
    move(Point(100, 0))
    arc(Point(0, 0), 100, 0, pi/4)
    arc(Point(0, 0), 200, pi/4, pi/2)
    strokepath()
end
@svg begin
    background("black")
    sethue("white")
    for j = 0:15
        line(100 * (-1)^j, -200 + j * 20)
    end
    strokepath()
end

move a line pracují v souřadnicích relativních k počátku v bodě (0,0)(0,0). K dispozici máme i rmove a rline, které pracují relativně k aktuální poloze.

@svg begin
    background("black")
    sethue("white")

    circle(Point(-100, 0), 5, action = :fill)
    circle(Point(100, 0), 5, action = :fill)
    
    move(Point(-100, 0))
    for j = 0:8
        rline(Point(50*cos(j*2pi/8), 50*sin(j*2pi/8)))
    end
    closepath()
    
    move(Point(100, 0))
    for j = 0:8
        rline(Point(50*cos(j*2pi/8), 50*sin(j*2pi/8)))
    end
    closepath()
    
    strokepath()
end

Makra se hodí k rychlému prozkoumávání a experimentování v Notebooku nebo editoru (VS Code umí obrázky správně zobrazovat).

Co když chceme výsledný obrázek exportovat?

function on_circle(j::Integer, n::Integer, R::Real)
    return Point(R*cos((j-1)*2pi/n), R*sin((j-1)*2pi/n))
end

function my_colorful_ngon(n::Integer, filename::AbstractString; R::Real = 200)
    n < 3 && error("At least 3 vertices are needed!")
    R <= 0 && error("R has to be positive!")

    Drawing(500, 500, filename)
    origin()

    for j = 1:n
        setcolor(rand(3)...)
        poly([ Point(0,0), on_circle(j-1, n, R), on_circle(j, n, R) ], :fill)
    end

    finish()
    preview()
end
my_colorful_ngon (generic function with 1 method)
my_colorful_ngon(3, "3.svg")
my_colorful_ngon(8, "8.svg")

Řešení některých příkladů

K Julia:

function julia(c::Complex{Float64}, R::Float64, reals, imags; imax::Int64=1000)
    img   = zeros(length(reals), length(imags))
    fc(z) = z^2 + c
    
    for j in axes(reals, 1), k in axes(imags, 1)
        n = 0
        z = reals[j] + 1im*imags[k]
        
        while n <= imax && abs(z) < R
            z = fc(z)
            n += 1
        end
        
        img[j, k] = n
    end
    
    return img
end