Hlavní navigace

Fraktály v počítačové grafice X

Pavel Tišnovský 28. 12. 2005

V desátém pokračování seriálu o fraktálech používaných (nejen) v počítačové grafice se budeme zabývat velmi zajímavým tématem - zobrazováním map vybraných dynamických systémů v komplexní rovině. Mezi známé dynamické systémy, se kterými se v komplexní rovině pracuje, patří i Juliovy množiny a k nim příslušná Mandelbrotova množina, jejíž obrázek je všeobecně považován (spolu se slavnou Newellovou čajovou konvicí) za jeden ze symbolů počítačové grafiky.

Obsah

1. Mapy dynamických systémů v komplexní rovině
2. Stručná historie Juliových množin
3. Komplexní parabola
4. Trojúhelníková nerovnost pro funkci komplexní paraboly
5. Definice Juliových množin
6. Vykreslení Juliových množin
7. Obsah dalšího pokračování tohoto seriálu

1. Mapy dynamických systémů v komplexní rovině

Pro studium dynamických systémů je mnohdy důležité prozkoumat jejich dynamické chování, to znamená postupnou změnu stavu systému v čase. V principu je možné stavy, v nichž se systém nacházel v určitém časovém intervalu, vizualizovat dvěma způsoby: animací nebo současným zobrazením stavového prostoru. Stav systému je popsán, jak již víme z předchozích pokračování tohoto seriálu, stavovým vektorem, je tedy nutné určitým způsobem vykreslit tento vektor. Pokud má stavový vektor jednu, dvě nebo tři složky, je vykreslení jednoduché – složky stavového vektoru určují pozici bodu v rovině či prostoru. Stav systému je potom znázorněn těmito body, které mohou být v případě potřeby pospojovány úsečkami, jak tomu bylo u Lorenzova a Rosslerova atraktoru. Je-li složek stavového vektoru více, je možné hodnotu zbývajících složek považovat například za barvu bodu či jeho průhlednost. Tento přímočarý způsob vykreslování jsme si již ukázali v předchozích částech tohoto seriálu na mnoha demonstračních příkladech. V následujícím textu si však ukážeme další možný způsob vizualizace dynamického systému, který je z hlediska počítačové grafiky velmi zajímavý a přínosný.

Fraktály 10 - 1

Obrázek 1: Dynamický systém v komplexní rovině

Místo drah orbitů budeme zobrazovat takzvané mapy dynamických systémů. Mapou je zde myšlen rastrový obrázek (pixmapa), ve kterém barva každého pixelu odpovídá vybranému stavu dynamického systému v bodě, jehož souřadnice jsou do pixelu mapovány – pixmapa je tedy „položena“ do komplexní roviny tak, že každému pixelu z pixmapy můžeme jednoznačně přiřadit jedinou komplexní hodnotu, tj. reálnou a imaginární složku. Zdaleka nejpoužívanějším způsobem vizualizace dynamických systémů v komplexní rovině je zvýraznění počtu iterací, tj. počtu opakování funkce dynamického systému do té doby, než je splněna nějaká předem známá podmínka. Naše povídání o dynamických systémech v komplexní rovině začneme popisem Juliových množin, z nichž si posléze (v dalších pokračováních tohoto seriálu) odvodíme pravidla pro tvorbuMandelbro­tovy množiny.

Fraktály 10 - 2

Obrázek 2: Další typ dynamického systému v komplexní rovině

2. Stručná historie Juliových množin

Juliovy množiny jsou pojmenovány po francouzském matematikovi, který se jmenoval Gaston Julia. Tento matematik se začal v roce 1917 (konkrétně to bylo v nemocnici, kde se za první světové války léčil z úrazu) zabývat problematikou analýzy chování funkce komplexní paraboly – viz další kapitola. Posléze na této problematice pracovali jak Gaston Julia, tak i další francouzský matematikPierre Fatou, podstatná část jejich společné vědecké práce je datována do roku 1920. Na tyto práce, které byly prakticky zapomenuty (převažoval totiž mylný názor, že do matematiky nic podstatného nepřináší), navázal v osmdesátých letech minulého století Benoit B. Mandelbrot, který problematiku analýzy komplexní paraboly spojil s fraktální geometrií (po Mandelbrotovi se tímto fenoménem zabýval zejména A. Douday a J. Hubbard). Po Gastonu Juliovi jsou pojmenovány Juliovy množiny a po Pierre Fatouovi jejich speciální podtyp, takzvaný Fatouův prach (Fatou dust). To však platí zejména pro evropskou literaturu, v americké se tento podtyp Juliových množin nazývá zobecněný Cantorův prach (Cantor dust). Nyní si však pojďme funkci komplexní paraboly popsat.

Fraktály 10 - 3

Obrázek 3: Pohled na Juliovu množinu

3. Komplexní parabola

Funkce komplexní paraboly vypadá následovně:

y=f(x)=x2+c

Gaston Julia se zajímal především o iterační proces, ve kterém je výsledek výpočtu v jednom kroku použit jako vstupní parametr v kroku následujícím, tj. funkce komplexní paraboly je zde „zapojena“ do zpětné vazby (proměnné x a y jsou nahrazeny řadou):

zn+1=f(zn)=z2+c

přičemž znc leží v komplexní rovině. Pokud by zn a c ležely v oboru reálných čísel, jednalo by se o funkci klasické paraboly, nejznámější je samozřejmě případ, kdy jec rovno nule. Pokud by se i v komplexní rovině zavedla nulová konstanta c, došli bychom k následující posloupnosti:

z0 → z02 → z04 → z08 → …

Pro tuto posloupnost je možné velmi jednoduše zjistit, jaký tvar bude mít atraktor dynamického systému tvořený právě funkcí komplexní paraboly:

  1. Pokud je velikost počáteční hodnoty z0 ostře menší než 1, tj. platí vztah |z0|<1, je atraktorem tohoto dynamického systému jediný bod – počátek komplexní roviny, tj. „komplexní nula“ 0+0i. Posloupnost hodnot zn v tomto případě směřuje k počátku, při zobrazení a propojení bodů v posloupnosti bychom viděli spirálu. Při každé iteraci se zdvojnásobí úhel φ komplexního čísla (to totiž může být vyjádřeno také jako |z|e), pro některé hodnoty z tedy nastává periodické opakování hodnoty úhluφ.
  2. Pokud je velikost počáteční hodnoty z0 ostře větší než 1, tj. platí vztah |z0|<1, je atraktorem nekonečno (komplexní rovina se pro tyto účely mapuje na takzvanou Riemannovu sféru (kouli), takže nekonečno je tvořeno jedním bodem koule – pólem; druhým pólem je komplexní nula). I v tomto případě tvoří posloupnost hodnotzn spirálu a samozřejmě i zde dochází při iterativním výpočtu k postupnému zdvojnásobování hodnoty úhlu φ.
  3. Pokud je velikost počáteční hodnoty z0 přesně jednotková, je atraktorem kružnice se středem v počátku komplexní roviny a poloměrem rovným jedné.

Jak je z předchozích tří bodů patrné, je tvar atraktoru v případě, že c je nulové, velmi jednoduše zjistitelný. Obtížnější situace nastane v případě, že je hodnota c nenulová. Zde se ukazuje, že pro některé hodnoty c a počáteční hodnoty z0 nenímožné zjistit, zda posloupnost zn bude konvergovat k nějakému fixnímu bodu, divergovat či zda půjde o cyklickou posloupnost. Navíc se ukazuje, že mapa bodů, které nedivergují, tvoří množinu s fraktálními vlastnostmi – soběpodobností a nezávislostí na změně měřítka.

Fraktály 10 - 4

Obrázek 4: Modifikovaná Juliova množina

4. Trojúhelníková nerovnost pro funkci komplexní paraboly

Při zjišťování, zda posloupnost zn+1=f(zn)=z2+c diverguje, můžeme využít následujícího tvrzení:

Pokud platí |c|<=|z| a současně |z|>2, pak posloupnost diverguje. Důkaz tohoto tvrzení je založen na trojúhelníkové nerovnosti – viz pátý ilustrační obrázek a sekvence odvození:

Fraktály 10 - 5

Obrázek 5: Trojúhelníková nerovnost pro funkci komplexní paraboly

|z| > 2 ^ |z2+c| >= |z2| – |c|     (musí platit pro každý trojúhelník)

|z| > 2 ^ |z2| > 2|z| → |z2 + c| > 2|z| – |c| (první důležitý mezivýsledek)

|z| > =|c| → |z| + |z| – |z| >= |c|
2|z| – |z| >= |c| → 2|z| – |c| >= |z| (druhý důležitý mezivýsledek)
|z| > 2 ^ |z| >= |c|

|z2+c| > 2|z| – |c| >= |z| (porovnání obou mezivýsledků)
|z2+c| > |z| (posloupnost opravdu diverguje)
 

Z předchozích vztahů tedy můžeme odvodit následující pravidla:

  1. Pokud je velikost zn větší jak 2, tj. |zn|>2, posloupnost začíná divergovat, a to bez ohledu na hodnoty dosažené v předchozích iteracích i hodnoty c (velikost|c| by však měla být menší než 2, viz další pravidlo). Laicky řečeno, takto velké zn „přebije“ všechny ostatní vlivy v iteračním výpočtu.
  2. Pokud je velikost konstanty c větší jak 2, tj. |c|>2, posloupnost opět diverguje, bez ohledu na hodnoty zn.
  3. Pokud je velikost konstanty c nulová, platí pro posloupnost vztahy uvedené v předchozí kapitole, tj. v závislosti na velikosti |z0| posloupnost konverguje k nule, diverguje nebo se nachází na jednotkové kružnici.
  4. V ostatních případech není možné analyticky dokázat, zda posloupnost konverguje, diverguje či osciluje, proto je nutné provádět další iterace. Těch však může být nekonečně mnoho, pro některé body nelze konvergenci či divergenci v konečném čase rozhodnout, což je velmi zajímavý výsledek, který souvisí s fraktálními charakteristikami Juliových množin.

5. Definice Juliových množin

Juliovy množiny jsou vytvářeny, podobně jako později popsaná Mandelbrotova množina a další typy fraktálů, pomocí jednoduchého dynamického systému, založeného na postupné iteraci výše zmíněné funkce komplexní paraboly:

zn+1=zn2+c

kde proměnné zn a c leží v komplexní rovině. Iterační proces začíná probíhat se startovní hodnotou z0, která v případě Juliových množin reprezentuje pozici bodu v komplexní rovině. Komplexní hodnota c je zvolena libovolně a pro všechny počítané body v jednom obrazci zůstává konstantní. Praktický význam však mají takové konstanty c, jejichž velikost nepřekročí hodnotu 2, což vyplývá ze závěrů napsaných v předchozí kapitole.

Juliovy množiny jsou definovány jako množiny všech komplexních čísel z0, pro které posloupnost zn nediverguje, tj. posloupnost buď konverguje nebo osciluje:

Fraktály 10 - 6

kde Pcn(0) značí hodnotu zn pro danou konstantu c a z0.

Fraktály 10 - 7

Obrázek 7: Další modifikace Juliovy množiny

6. Vykreslení Juliových množin

Vykreslení Juliových množin je poměrně přímočaré, jak je tomu ostatně i u dalších typů fraktálních útvarů. Dvě vnější programové smyčky zajišťují generování počátečních souřadnic bodů z0, které jsou posléze iterovány v nejvnitřnější smyčce. Každá iterace probíhá tak, že se postupně počítají hodnoty z1, z2 atd., do chvíle, než velikost nějaké hodnoty zn překročí hodnotu 2, protože o této hodnotě již z výše uvedeného důkazu víme, že představuje mez pro divergenci. Problém nastává v případě, že divergenci není možné zjistit v konečném počtu kroků. Zde nastupuje malý trik v tom, že vnitřní iterační smyčka má nastaven maximální počet opakování. Pokud se ani po tomto počtu opakování nezjistí divergence, je posloupnost považována za konvergující, což je samozřejmě z matematického hlediska nesmysl (protože hned v dalším kroku se například může zjistit divergence). Podle počtu provedených iterací je obarven pixel, jehož souřadnice odpovídají souřadnici počátečního bodu z0. Počátečními parametry pro tento výpočet je maximální počet iterací, a také hodnota komplexní konstanty c, která je v následující proceduře uložena do proměnných cx (reálná složka) a cy (imaginární složka).

//-----------------------------------------------------------------------------
// Procedura provádějící překreslení Juliovy množiny se standardním výpočtem
//-----------------------------------------------------------------------------
void recalcJulia( pixmap *pix,                  // pixmapa pro vykreslování
                  int    maxiter)               // maximální počet iterací
{
    double zx, zy, zxx;                         // složky komplexní proměnné Z
    double zx0, zy0;
    double cx=0.0, cy=0.9;                      // složky komplexní konstanty C
    int    x, y;                                // počitadla sloupců a řádků v pixmapě
    int    iter;                                // počitadlo iterací

    zy0=MINY;
    for (y=0; y<pix->height; y++) {             // pro všechny řádky v pixmapě
        zx0=MINX;
        for (x=0; x<pix->width; x++) {          // pro všechny pixely na řádku
            zx=zx0;                             // nastavit počáteční hodnotu Z(0)
            zy=zy0;
            for (iter=0; iter<maxiter; iter++) {// iterační smyčka
                if (zx*zx+zy*zy>4.0) break;     // kontrola překročení meze divergence
                zxx=zx*zx-zy*zy+cx;             // výpočet Z(n+1)=Z(n)^2+C
                zy=2.0*zx*zy+cy;
                zx=zxx;
            }
            putpixel(pix, x, y, iter, iter, iter); // vykreslení pixelu
            zx0+=(MAXX-MINX)/pix->width;        // posun na další bod na řádku
        }
        zy0+=(MAXY-MINY)/pix->height;           // posun na další řádek
    }
} 

Vzhledem k tomu, že je generování Juliových množin obecně pomalé (zejména při nastavení vyššího počtu iterací), je zapotřebí provést optimalizaci vnitřní iterační smyčky – alespoň na úrovni céčkovského kódu, i když mnohé programy obsahují vnitřní iterační smyčky optimalizované ručně v assembleru. Při pohledu na první verzi iterační smyčky je patrné, že se používá šest operací násobení a čtyři operace sčítání a odečítání. Vhodnou úpravou výrazů a použitím dvou pomocných proměnných je možné snížit počet operací násobení na pouhé čtyři, jak je to ukázáno na upravené proceduře zobrazené níže. V této proceduře se používají pomocné proměnné nazvanézx2 a zy2, které obsahují druhou mocninu reálné (zx) a imaginární (zy) složky proměnné z. Současně jsou prohozeny výrazy pro výpočet nové hodnoty zx a zy, čímž ušetříme použití pomocné proměnné zxx (nová hodnota proměnné zx totiž není po úpravách závislá na staré hodnotě zy).

//-----------------------------------------------------------------------------
// Procedura provádějící překreslení Juliovy množiny se zrychleným výpočtem
//-----------------------------------------------------------------------------
void recalcJulia( pixmap *pix,                  // pixmapa pro vykreslování
                  int    maxiter)               // maximální počet iterací
{
    double zx, zy, zx2, zy2;                    // složky komplexní proměnné Z a Z^2
    double zx0, zy0;
    double cx=0.0, cy=0.9;                      // složky komplexní konstanty C
    int    x, y;                                // počitadla sloupců a řádků v pixmapě
    int    iter;                                // počitadlo iterací

    zy0=MINY;
    for (y=0; y<pix->height; y++) {             // pro všechny řádky v pixmapě
        zx0=MINX;
        for (x=0; x<pix->width; x++) {          // pro všechny pixely na řádku
            zx=zx0;                             // nastavit počáteční hodnotu Z(0)
            zy=zy0;
            for (iter=0; iter<maxiter; iter++) {// iterační smyčka
                zx2=zx*zx;                      // zkrácený výpočet druhé mocniny složek Z
                zy2=zy*zy;
                if (zx2+zy2>4.0) break;         // kontrola překročení meze divergence
                zy=2.0*zx*zy+cy;                // výpočet Z(n+1)
                zx=zx2-zy2+cx;
            }
            putpixel(pix, x, y, iter, iter, iter); // vykreslení pixelu
            zx0+=(MAXX-MINX)/pix->width;        // posun na další bod na řádku
        }
        zy0+=(MAXY-MINY)/pix->height;           // posun na další řádek
    }
} 
Fraktály 10 - 8

Obrázek 8: Ještě jedna modifikace Juliovy množiny

7. Obsah dalšího pokračování tohoto seriálu

V dalším pokračování tohoto seriálu si řekneme, do jakých kategorií se Juliovy množiny dělí. Také si ukážeme praktické příklady vykreslování Juliových množin na několika demonstračních příkladech.

Našli jste v článku chybu?

3. 1. 2006 9:10

Je pravda, ze ja jsem se o architekturu modernich procesoru prestal zajimat uz od Pentia - o nem jsem jeste prelouskal par clanku venovanych optimalizacim, ale potom jsem si uvedomil, ze nebudu mit cas a chut provadet rucni optimalizaci pro jednotlive pipeliny (je to taky zbytecne, protoze kazda nova generace procesoru je v tomto odlisna). Od te doby jsem prakticky na assembler nesahl :-|

Pravda je, ze pokud prekladac cecka preklada podle normy, mel by vysledek kazdeho vyrazu prevest na dany da…

3. 1. 2006 9:01

Taky me to po me predchozi odpovedi napadlo, protoze jsem si az potom vsiml, ze v Rootovske strance chybi tagy pro konce radku :-). S timto vsak nemuzu nic udelat, muj puvodni XHTML kod (s "pre" a samozrejme pouze s CR) je automaticky prevaden redakcnim systemem, ktery pouziva "code" a samozrejme nejake CSS k tomu, aby to bylo podsviceno bile a pripadne se zobrazovaly scrollbary. Asi zbyva bud pockat do zitrka na dalsi pokracovani :-) nebo ve zdrojaku XHTML, ktery Root posila…
Podnikatel.cz: Přehledná titulka, průvodci, responzivita

Přehledná titulka, průvodci, responzivita

Lupa.cz: Proč firmy málo chrání data? Chovají se logicky

Proč firmy málo chrání data? Chovají se logicky

Podnikatel.cz: EET: Totálně nezvládli metodologii projektu

EET: Totálně nezvládli metodologii projektu

Podnikatel.cz: Vrátí zvýhodnění, ale výrazně omezí paušály

Vrátí zvýhodnění, ale výrazně omezí paušály

Lupa.cz: Google měl výpadek, nejel Gmail ani YouTube

Google měl výpadek, nejel Gmail ani YouTube

Lupa.cz: Co se dá měřit přes Internet věcí

Co se dá měřit přes Internet věcí

120na80.cz: Rakovina oka. Jak ji poznáte?

Rakovina oka. Jak ji poznáte?

Měšec.cz: Finančním poradcům hrozí vracení provizí

Finančním poradcům hrozí vracení provizí

Podnikatel.cz: Babiš: E-shopy z EET možná vyjmeme

Babiš: E-shopy z EET možná vyjmeme

Měšec.cz: Air Bank zruší TOP3 garanci a zdražuje kurzy

Air Bank zruší TOP3 garanci a zdražuje kurzy

Podnikatel.cz: Babiše přesvědčila 89letá podnikatelka?!

Babiše přesvědčila 89letá podnikatelka?!

Měšec.cz: Kdy vám stát dá na stěhování 50 000 Kč?

Kdy vám stát dá na stěhování 50 000 Kč?

Lupa.cz: Teletext je „internetem hipsterů“

Teletext je „internetem hipsterů“

DigiZone.cz: Rádio Šlágr má licenci pro digi vysílání

Rádio Šlágr má licenci pro digi vysílání

Vitalia.cz: Jsou čajové sáčky toxické?

Jsou čajové sáčky toxické?

Měšec.cz: Jak vymáhat výživné zadarmo?

Jak vymáhat výživné zadarmo?

Podnikatel.cz: Podnikatelům dorazí varování od BSA

Podnikatelům dorazí varování od BSA

Root.cz: Vypadl Google a rozbilo se toho hodně

Vypadl Google a rozbilo se toho hodně

Vitalia.cz: To není kašel! Správná diagnóza zachrání život

To není kašel! Správná diagnóza zachrání život

Lupa.cz: Insolvenční řízení kvůli cookies? Vítejte v ČR

Insolvenční řízení kvůli cookies? Vítejte v ČR