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 GRZde 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 GRTvorbu 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
xayreprezentující souřadnice bodů, - vektor
xa funkce pro výpočet -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ě metodatitle),xlabel: titulek grafu (alternativně metodaxlabel),ylabel: titulek grafu (alternativně metodaylabel),xlim,ylim: rozsah jednotlivých os (alternativně metodyxlimaylim),
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 (vizGR.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í: ???
- 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))xyscfigure()
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)xyscplot3
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
endfigure(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)zsheatmap
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í zadané předpisem , pro pevně zadanou konstantu . Zkoumáme pak množinu bodů , pro které je posloupnost 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é a maximální počet iterací a pro z nějaké zvolené množiny (typicky mřížka) naměřme, kolik iterací musíme provést, než se dostaneme dále než od počátku.
Proveďte tento výpočet pro a .
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)
endres = 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 PyPlotres = 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 PyPlotPro 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 PyPlotplot
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 , , pomocí maticové exponenciály. Řešením byla vektorová funkce .
Konkrétně jsme zvolili
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)
# ...
endts = 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 PGFPlotsXusing PGFPlotsXTakř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 CairoMakiexs = 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
enddata = my_prng(10^6, 0.7500000001);Jak tato posloupnost vypadá?
using CairoMakiefig, 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, CairoMakiexs = 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))
end2. 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, Colors2.1 Vytvoření obrázku
img_data = rand(400, 600)
img_data[1:50, 1:50] .= 1
img_dataimg = 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 FileIOsave("test.jpg", img)Nahrávání obrázků
using FileIO, ImageMagick, ImageIOimg_Gauss = load("../homeworks/B211/files/gauss_portret.jpg")typeof(img_Gauss)dump(img_Gauss)Všimněte si speciálního datového typu:
?N0f8Jeden pixel je barevně reprezentován trojicí čísel (RGB).
p = img_Gauss[1, 1]typeof(p)p.rp.gp.bimg_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()
endPozor, 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()
endmove a line pracují v souřadnicích relativních k počátku v bodě .
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()
endMakra 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()
endmy_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