Nem lehet szó nélkül elmenni a 20. század hivatalosan is megválasztott (rangsorban nyolcadik) legjelentősebb algoritmusa mellett, ami gyors Fourier transzformáció névre hallgat (kiejtés: furié). Valamilyen oknál fogva ez kimaradt az egyetemi éveimből, holott csak a felhasználásait tekintve is forradalminak számít. Jelfeldolgozástól kriptográfiáig rengeteg területen használják; természetesen ezen viszonylag széles skálából most a számítógépes grafikában való alkalmazásairól fogok írni, konkrétan a valósidejű óceán megjelenítésről.

A címben szándékosan nem "szimulációt" írtam, mert bár a módszer egy mérhető adatokon alapuló szimuláció, annak csak a Fourier-analízis felőli részével szeretnék most foglalkozni, nem pedig a sokkal nehezebb, fizikai interakciókat is leíró Navier-Stokes egyenletekkel (majd talán egy következő cikkben). Az irodalmak közül kiemelném Jerry Tessendorf nevezetes dolgozatát, a használt összehasonlítási alapként pedig az nVidia implementációját.

A cikk matematikai része első olvasatra viszonylag nehéz, és engem sem segít az, hogy a legtöbb forrás a jelfeldolgozás oldaláról írja le a dolgokat, mellette nem kis kavarásokat okozva. A későbbi bekezdések a már kész algoritmus felhasználásáról szólnak egy minél realisztikusabb óceánfelület megjelenítéséhez. Ez utóbbiak többnyire időigényes programozási feladatok, ezért csak mértékletesen foglalkozok velük.

Aki úgy gondolja, hogy eleget tud a komplex számokról, az átugorhatja ezt a bekezdést. Történelemleckeként annyit, hogy már az ókori görögöknél (Hérón) is felmerültek olyan problémák, ahol negatív számok gyökét kellett volna meghatározni, de mint ahogy az irracionális számok koncepcióját, úgy ezt is kihajították a...Musaeum ablakán...

A 16. században már nagyobb fejtörést okozott, hogy az olyan harmadfokú egyenletek megoldóképletében, amelyeknek három valós gyöke van, megjelennek olyan részkifejezések, amik negatív számok valamilyen gyökét (radical) tartalmazzák. A probléma ugyan megkerülhető trigonometrikus függvények használatával, de az igényességre való törekvés végül elvezetett a komplex számok kidolgozásához Rafael Bombelli által.

(megj.: Bombelli valójában csak annyit mutatott meg, hogy hogyan lehet √-1-el számolni. A komplex számok itt bemutatott modern definícióját Hamilton dolgozta ki, illetve tovább általánosítva bevezette a kvaterniók fogalmát: ℍ = ℂ × ℂ)

A komplex számok halmazát matematikailag úgy vezetjük be, hogy vesszük a valós számok halmazának önmagával vett direkszorzatát (hasonlóan, mint a 2D vektorok esetében). Ez jelölésben ℂ = ℝ × ℝ. A vektoroktól való különbség a felírásban, illetve az elvégezhető műveletekben jelentkezik. Utóbbiak közül rögtön a szorzást emelném ki, egyelőre a rendezett pároknál megszokott jelöléssel:

Ezzel a definícióval élve, ha vesszük a (0, 1) komplex szám önmagával vett szorzatát (négyzetét), akkor meglepő módon (-1)-et kapunk. Jelüljük az említett komplex számot i-vel (fizikusok j-vel), ekkor egy tetszőleges komplex szám


alakban írható, ahol a és b valós számok, és a-t a komplex szám valós részének, b-t pedig a képzetes részének nevezzük. A jobb oldali jelölés előnye, hogy a valós számokból ismert képleteinket ugyanúgy alkalmazhatjuk, nem elfelejtve, hogy i² = -1. A továbbiakban ezt a fajta jelölést fogom használni (bár megemlíteném, hogy papíron az előbbi időnként kényelmesebb).

Amiért ez megtehető, az azzal indokolt, hogy a komplex számok halmaza is ún. (szám)testet (field) alkot a fent bemutatott szorzással, illetve a tagonként vett összeadással. A műveletek kommutativitása, asszociativitása, disztributivitása illetve invertálhatósága némi belegondolással levezethető. Egy olyan művelet van, amit valós számoknál nem használunk, ez a konjugált:


A konjugáltat *-al is szokták jelölni; későbbi bekezdésekben én is ezt fogom használni (tipográfiai szempontból átláthatóbb).

Jogosan merül fel a kérdés, hogy vajon a komplex számok halmaza izomorf-e a 2D vektorokkal? Ameddig halmazként tekintünk rájuk, addig igen, a műveleteket figyelembe véve viszont értelemszerűen nem. Természetesen ez nem akadályoz meg senkit abban, hogy komplex számokat helyvektorokként vizualizáljon, különösen amikor az ún. n-edik komplex egységgyökökről van szó (zⁿ = 1). Ezeknek a (diszkrét) Fourier transzformációval való kapcsolatáról majd később írok.

Első példaként határozzuk meg a 4x² + 1 = 0 egyenlet gyökeit. Az általános iskolában tanult másodfokú egyenlet megoldóképlete alapján a diszkrimináns -16, amire akkor azt mondtuk, hogy az egyenletnek nincs (valós) megoldása. Felvértezve azonban a komplex számokkal √(-16)-ből kiemelhetünk i² = -1-et, így ±i/2-t kapva egyszerűsítés után.

Ezek után megfogalmazható az algebra alaptétele, miszerint minden n-edfokú polinomnak pontosan n darab gyöke van (multiplicitást is beleértve).

A témához kapcsolódóan két probléma foglalkoztatta a 18. század vezető matematikusait: a rezgő húr probléma (d'Alembert, Euler, Bernoulli), illetve a hővezetés egyenlete (Fourier).

A rezgő húr problémát konyhanyelven úgy lehet megfogalmazni, hogy két pont között kifeszítünk egy rugalmas húrt és megpendítjük. A kérdés az, hogy hogyan lehet matematikailag leírni a húr mozgását (rezgését). Hosszas vita után Bernoulli a kortárs matematikusok eredményeit felhasználva arra a következtetésre jutott, hogy a rezgést súlyozott szinusz függvények összege (mint végtelen sorösszeg) írja le. Ezt a megoldást akkor még nem fogadták el széles körben (többek között azért, mert tisztázatlan volt pl. a függvény és az integrálás fogalma).

Több mint 50 évvel később Joseph Fourier a hővezetés egyenletét vizsgálva ugyanarra az eredményre jutott, mint annak idején Bernoulli. Szintén konyhanyelven: adott egy a környezetétől elszigetelt homogén rúd, a feladat pedig a rúd belsejében történő hőmérsékletváltozások matematikai leírása. Az eredményei alapján egy meglepő állítással állt elő: minden 2π szerint periodikus függvény felírható súlyozott szinusz és koszinusz függvények végtelen sorösszegeként. Ezt Lagrange visszautasította azzal, hogy "a sarkokat nem lehet szinuszokkal előállítani".

A mondjuk így "nagy előrelépés" hosszúnevű Dirichlet 1829-ben írt nevezetes dolgozatához fűződik, amiben egyrészt letisztázta a függvény, mint fogalom pontos meghatározását, illetve immár egzakt matematikai felírással és konkrét bizonyításokkal támasztotta alá a Fourier sorok elméletét.

Amire egyébként több ekvivalens definíció is létezik; én abból fogok kiindulni, amit az egyetemen is tanultam (bár ennél lényegesen bonyolultabb formában). Vegyük az f: ℝ → ℝ valamilyen P ∈ ℝ szerint periodikus függvényt, ami integrálható az [x₀, x₀ + P] intervallumon valamilyen x₀ ∈ ℝ-re. Ekkor az f függvény Fourier részletösszeg-függvényének (röf, igen...) nevezzük az

összeget, ahol



(megj.: s_N(x) nem a legszerencsésebb jelölés, de ennek megindoklása csak tovább bonyolítaná a cikket)

Ha N helyére végtelent írunk, akkor az így kapott kifejezés az f függvény Fourier sora. Felmerül a kérdés, hogy a sor milyen feltételek mellett konvergens (azaz N → ∞ esetén mikor és hogyan állítja elő a függvényt).

A válasz az, hogy "attól függ milyen függvény", illetve "attól függ milyen konvergencia". Ez a harmonikus analízis témaköre, ami elég körültekintően kivesézi az egyes eseteket, de ehhez a cikkhez csak néhány egyszerűbbet említek meg:

• ha a függvény folytonos^(*) és deriválható az x pontban, akkor a sor abban a pontban konvergens
• ha a függvény valamelyik L^p térből való (1 < p < ∞), akkor a sor majdnem mindenhol konvergens a p-normában

Utóbbihoz kapcsolódik a Riesz-Fisher tétel, ami kimondja, hogy ha a függvény L²-beli (négyzetesen integrálható), akkor a Fourier sora konvergens a 2-normában (az eltérés négyzetének integrálja tart nullához). Vigyázat: az, hogy a sor konvergens még nem jelenti, hogy elő is állítja a függvényt! ^(*)A lenti ábrán is látható lesz, hogy a törési pontokban a sor túlleng a függvényen (ún. Gibbs jelenség), és ez sajnos N = ∞ esetén sem tűnik el. Ezeken a helyeken a sor a függvény jobb- és baloldali határértékének átlagához konvergál, amennyiben létezik jobb- és baloldali deriváltja.

Más okból is célszerű legalább L²-beli függvényeket választani, ugyanis az előbb említett feltétel p = 1 esetén nem igaz (értsd: nem elég, hogy a függvény "csak" integrálható). Konstruálható például olyan tisztán L¹-beli függvény, aminek a Fourier sora divergál.

De inkább nézzük meg egy példán keresztül, hogy mire lehet számítani. Vegyük mondjuk a fűrészfog függvényt, amit az egyszerű f(x) = x - floor(x) képlettel lehet leírni (egyébként meg pl. szintetizátorok használják hangelőállításra). A függvény P = 1 szerint periodikus, illetve a [0, 1] intervallumon integrálható, itt ugyanis elég úgy kezelni, mintha a g(x) = x függvény lenne. Az eredményt N = 4 és N = 16-ra mutatom:


Ilyet már láttunk, amikor például convolution shadow map-et simítottunk el a Fourier sora segítségével. Mivel a fűrészfog függvény páratlan (f(-x) = -f(x)), az a_k együtthatók mind nullák, a békák b_k-k pedig meghatározhatóak parciális integrálással. Magát a sort nem képlettel írom fel, hanem gnuplot kóddal (hogy ilyet is lássatok):


Mint mondtam, shadow mapping-ben már láttuk egy gyakorlati alkalmazását a Fourier soroknak, de jogos a kérdés, hogy általános esetben miért jó egy ehhez hasonló előállítás (már ha az nem elég, hogy a rezgő húrt illetve a hővezetés egyenletét leírja). Mivel szinusz/koszinusz függvények összegéről van szó, a sor kapásból teljesít egy rakatnyi hasznos tulajdonságot (pl. végtelenszer folytonosan differenciálható), amik által az eredeti függvény bizonyos viselkedéseit/tulajdonságait is meg lehet jósolni (és ez a méréselméletben különösen hasznos, ami által eljutunk majd az óceánfelület analitikus modelljeihez).

Utolsó megjegyzésként, ha a fenti Fourier röf definíciót gondolatban szétbontjuk, akkor az vehető észre, hogy az egyes tagok frekvenciái az f₀ = 1/P alapfrekvencia, pontosabban az ω₀ = 2π/P alap körfrekvencia egész számú többszörösei. A szumma így előálló k-adik tagját nevezzük k-adik felharmonikusnak.

Az egyik részcél az lesz, hogy (megint konyhanyelven) a "sütiből előállítsuk a receptet". Ez matematikailag azt jelenti, hogy valamilyen (analóg vagy digitális) jelből meghatározzuk az egyes felharmonikusok amplitúdóit és fázisait a frekvencia függvényében. Az így kapott diagramo(ka)t nevezzük spektrumnak.

(megj.: egyéb irodalmak esetleg máshogy értelmezik a spektrum fogalmát, illetve korábban említettem, hogy egy mátrix sajátértékeit is szokás spektrumnak hívni; későbbi bekezdésekben majd kiderül, hogy van-e kapcsolat a kettő között)

A fenti felírás egyik problémája, hogy két (látszólag) különböző Fourier együttható van, amiket nem is mindig trivi kiszámolni. Célszerű lenne egy olyan alakban felírni a Fourier sort, hogy csak egy darab együttható legyen, ami cserébe viszont komplex. A fázist is szemléletesen lehet vizualizálni egy ilyen alakban, mint a komplex együttható és a valós számtengely által bezárt szöget. Az átíráshoz induljunk ki az Euler formulából:

Ha megnézzük a jobb oldalt, az a korábban már látott felírása egy komplex számnak, ami alapján a koszinusz felírható az e^ix valós részeként, a szinusz pedig a képzetes részeként. Ha hasonló módon felírjuk az Euler formulát e^-ix-re is, akkor a kapott két egyenletet összeadva illetve kivonva kapjuk, hogy:


Mielőtt ész nélkül behelyettesítenénk ezeket a fenti Fourier sor definícióba, érdemes alkalmazni a már említett körfrekvencia képletet (ω = 2π/P), amivel kicsit érthetőbben lehet felírni az átalakításokat:


Höfö ellenőrizni, hogy nem számoltam-e el (tudom, hogy ez is röf, de nem röfögök most már többet; ez nem elszámolás hanem elírás). A c_k együtthatók a következőképpen írhatóak fel kompaktabb alakban:


Remélem senki sem siklott át azon tény felett, hogy (mint említettem) a konjugáltat *-al jelölöm, másrészt, hogy az új szumma -N-től N-ig megy. Emiatt megjelentek negatív körfrekvenciák is a képletben, ami nem egy létező jelenség, hanem ennek a felírásnak egy tulajdonsága. Mivel a negatív/pozitív körfrekvenciákhoz tartozó értékek egymás konjugáltjai, összeadva őket az eredmény valós lesz.

A további magyarázkodások elkerülése végett itt írom le, hogy egy adott felharmonikus amplitúdója megkapható a c_k komplex szám hosszaként, a fázisszög pedig a képzetes és valós rész hányadosának arkusz tangenseként.

Egy pillanatra felejtsük el a "diszkrét" szót, és nézzük meg, hogy konkrétan hogyan néznek ki a felharmonikusok, illetve az (amplitudó)spektrum. Példának továbbra is jó lesz a fűrészfog függvény Fourier sora.

A bal oldali ábrán látható pirossal az eredeti Fourier sor (most mint jel), illetve annak felharmonikusai (amiket összeadva előáll). Amikor arra vagyunk kíváncsiak, hogy egy időtartománybeli jel tartalmaz-e egy adott frekvenciát, azt a frekvenciatartományba való transzformálásával tudjuk ellenőrizni. A jobb oldali ábra ezt a frekvenciatartománybeli felírást mutatja, mint az egyes felharmonikusok amplitúdóit.

(megj.: itt most az f₀ = 1 alapfrekvencia alapján rajzoltam meg a diagramot, de nyilván körfrekvenciára is teljesen hasonló)

A frekvenciatartományba való átalakítást írja le a Fourier transzformáció analóg (folytonos) jelekre, illetve a diszkrét Fourier transzformáció (DFT) digitális (mintavételezett) jelekre. Tegyük fel, hogy felvettünk valamilyen zenét, viszont valami interferencia miatt egy magas frekvenciájú sípolás került a felvételbe. Ilyenkor a spektrum előállításával ez megtalálható, eltüntethető, és egy inverz Fourier transzformációval megkapjuk a javított zenét.

A továbbiakban a DFT-vel fogok foglalkozni, mert a gyakorlatban ez dominál (a gép ezzel tud dolgozni). A definíció valamilyen {x[0]...x[N-1]} minták sorozatára (k ∈ [0, N - 1]):


A DFT tehát egy {x[k]} valós vagy komplex sorozatot alakít át az ugyannyi elemű {X[k]} komplex sorozatba. A műveletigényére ebből rögtön meg is mondható, hogy négyzetes (Ο(N²)). A másik amit észreveszünk, hogy felírható mátrix-vektor szorzásként is. Milyen transzformációkat is írtunk le mátrixokkal? Például lineárisakat (és valóban, ez egy lineáris transzformáció). Emellett invertálható is:


Különféle irodalmak különféle helyekre teszik az N-el való osztást; ennek megvan a magyarázata, de a fenti példából kiindulva világos lesz, hogy most a DFT-be kell az osztás és nem az inverzébe. Egyelőre vegyünk egy 2-periódusú szinusz hullámot. Ha jól számoltuk ki a DFT-t (N = 32 mintára), akkor a 0.5 alapfrekvenciában kell megkapjuk a hullám amplitúdóját (ugyanúgy a megfelelő X[k]-k hosszát véve):


(megj.: az említett negatív körfrekvenciák miatt kétszer kaptuk meg. A DFT-nek tehát csak az első N/2 tagja hordoz tényleges információt)

Egy kicsit arról, hogy hogyan is készítettem el az ábrákat. A bal oldal azt mutatja, hogy 4 másodpercig mintavételeztem, összesen 32 mintát, azaz a mintavételezési frekvencia 8 minta/másodperc. Innen a mintavételi időtartam h = 1/8 másodperc. A frekvenciatartomány lépésközét ebből lehet meghatározni a Δf = 1 / (N*h) képlettel.

Az is feltűnő, hogy a spektrum az említett helyeken kívül is nemnulla; ez a Fourier transzformáció egy "tulajdonsága", hogy gyakran Dirac delta függvényekkel fejezzük ki az eredményt (konyhanyelven "az ott valójában egy függőleges vonal, gyökér").

Néhány mondat a mintavételezésről: világos, hogy ha magas frekvenciákat tartalmazó jelet nem elég sűrűn mintavételezünk, akkor teljesen más eredményt kapunk, mint kellene (aliasing). A Nyquist mintavételezési tétel azt mondja, hogy a mintavételezési frekvencia nagyobb kell legyen a jel legmagasabb frekvenciájának kétszeresénél (Nyquist rate). A mintavételezési frekvencia fele (Nyquist limit) azt a maximális frekvenciát jelenti, amit még reprezentálni lehet. Például a CD lemeznek 44 kHz a mintavételezési frekvenciája, azaz legfeljebb 22 kHz-es felhangokig tárolja a felvett zenét.

(megj.: nem ez a baj vele, mert ennél magasabb frekvenciákat úgysem hallunk, hanem az, hogy az egyes minták 16 biten vannak tárolva, emiatt az elég nagy amplitúdókat levágja. Ha ezt kompenzálják, akkor pedig az apró részletek nem lesznek hallhatóak. Ez pontosan ugyanaz a probléma, mint grafikában az LDR vs. HDR)

Mint mondtam a DFT Ο(N²) komplex szorzás-összeadás párossal számolható ki, ami elég lassú (amennyiben részletes óceánfelületet szeretnénk majd megjeleníteni). Ordít viszont róla, hogy egy compute shader segítségével közel lineáris időre hozható:


A 60-as évek matematikusai kezüket-lábukat adták volna egy ilyen megoldásért, de mivel akkor még a GPU sem létezett, kénytelenek voltak más módszerekkel gyorsítani.

Kissé meglepő, hogy ezt az algoritmust Gauss is feltalálta a 19. század elején, de több szerencsétlen ok miatt nem terjedt el széles körben, ezért elfelejtődött. Az algoritmus újrafeltalálása James Cooley és John Tukey 1965-ös dolgozatában történt meg, így nem túl meglepő módon az egyik megoldást Cooley-Tukey FFT algoritmusnak hívják (de még ezen belül is van több változata).

Én most a radix-2 módszert fogom megmutatni (az nVidia implementációja a radix-8 módszert használja), ami azt jelenti, hogy kétpontos DFT-kből indul ki, illetve a műveletigénye Ο(N log₂N). Először nézzük meg, hogy mit is csinál egy kétpontos DFT. A képlet alapján könnyen rá lehet jönni, hogy (mivel most k, n ∈ [0, 1])

Ami alapján a kétpontos DFT:


Ezt hívják pillangó műveletnek és az alábbi folyamatábrával szokás vizualizálni (nem írtam el):


Egy általános pillangó művelet valamilyen a, b, α ∈ ℂ értékekre az a + αb és az a - αb értékeket adja meg. Az ábrán α = 1, ahogy az előbb kiszámoltam (egyébként majd az ún. twiddle factor fog odakerülni), az alul levő -1 pedig szintén az előbb kiszámolt e^-iπ (ez mindig -1 lesz).

Lépjünk tovább a négypontos DFT-re (és ennél tovább nem is fogok menni, mert nehéz megrajzolni a pillangó diagramokat). Hasonlóan mint előbb, kihasználva, hogy N = 4:


Amit kifejtve kapjuk, hogy:


És átrendezésekkel alakítsuk át pillangó műveletekké:


Ha összeszámoljuk ez összesen 4 darab pillangó művelet, tehát 8 darab komplex szorzás és összeadás (16 helyett). Az FFT ezek szerint úgy gyorsítja az algoritmust, hogy egyszerű manipulációkkal újrafelhasználja a már kiszámolt részeredményeket. Lássuk ennek a számolásnak a diagramját (figyeljétek bal oldalt az indexeket!):


Oké, szóval az átrendezések miatt a kiinduló sorozatot is át kell rendezni, de ha jobban megnézzük, ez az indexek bitjeinek megfordításával megtehető (és erre van beépített függvény GLSL-ben). Ezen kívül az ábrán most már a twiddle factor-ok általános jelölését használtam:


Amik egyébként az N-edik komplex egységgyökök és ráadásul N szerint periodikusak, így előre ki lehet őket számolni (emiatt szoktak nagyobb radixokkal dolgozni). Egy kicsit matematikaibb oldalról megközelítve, az FFT szétbontja a DFT-t a páros illetve páratlan indexű tagjaira, ami általánosan felírva (a páratlan szummából kiemelve a twiddle factor-t):


Az említett periodicitás miatt pedig:


A kettő együtt láthatóan egy pillangó művelet. Ha ugyanezt a szétbontást rekurzívan eljátsszuk a bal és jobb oldali szummákra, akkor log₂N lépés után már a kétpontos DFT-t kapjuk. A legintuitívabb implementáció tehát a rekurzió.

De mint tudjuk, shader-ekben nincs olyan, szóval egy iteratív megoldást kellene találni. Először megmutatom a wikipédián is megtalálható implementációt, felkommentezve (a bejövő adatot eleve az X[k] tömbbe rendeztem át, így csak azt fogom módosítani):

(megj.: pillangó csoport alatt az egymást keresztező pillangó műveletek halmazát értem)


(megj.: a jelölés ne zavarjon össze senkit; ez értelemszerűen nem ugyanaz az m és k, mint a képletben!)

Igen, ebben az algoritmusban vannak olyan dolgok, amiket nem magyaráztam el, például hogy az egyes szakaszokban lehet használni kisebb komplex egységgyököket is (számoljatok utána, ugyanaz fog kijönni). Emellett helyben végzi el az FFT-t, mint mondtam azzal az előfeltétellel, hogy a tömb már át van rendezve.

A valódi ok, amiért ezt az algoritmust megmutattam az, hogy némi módosítással triviálisan átírható compute shader-be. Ehhez a két belső ciklust (ami ha kiszámoljátok szakaszonként N/2 pillangó műveletet végez el) kellene egyesíteni egy darab ciklusba, ami viszont N műveletet végez el (amivel mellesleg nyilvánvalóbb is az Ο(N log₂N) műveletigény). Ez persze némi redundanciát jelent majd, de bőven meg fogja érni. Gondoljuk akkor végig, hogy hogyan lehetne ezt megcsinálni.

Az első logikus dolog, hogy a twiddle factor-t ne ciklussal számoljuk ki, hanem egy lépésben (és akkor a diagramra/képletre is jobban fog hasonlítani). A fenti algoritmus, mint mondtam egy adott szakaszra csak egyszer végzi el ugyanazt a pillangó műveletet; ez most meg fog változni annyira, ahányszor szerepel az eredeti felírásban.

Ki kell számolni tehát, hogy egy adott szakaszban (redundanciát is figyelembe véve) hol kezdődnek a pillangó műveletek. Ha ez megvan, akkor az előző szakasz eredményéből lehet kiszámolni az aktuális értékeket (azaz kell plusz egy tömb, a számolásokat pedig az eredeti és eközött pingpongoztatom el a jobb oldalig). Nézzük meg a kódot, érthetőbb:


Egy kis magyarázat: ha tudjuk, hogy a maradékos osztó kettőhatvány, akkor pozitív számokra a % m == a & (m - 1); ebből némi számolással levezethető, hogy miért a megfelelő értékeket kapjuk az indexekre. Mivel az FFT kettőhatvány méretű sorozatokon dolgozik, az eredmény garantáltan az X tömbben lesz. Ha esetleg az elemek száma nem kettőhatvány, akkor ki kell egészíteni az eredeti tömböt nullákkal.

Ezt a módosított algoritmust már könnyedén át lehet írni compute shader-be, amiben a belső ciklus eltűnik, a ciklusváltozót (l) pedig megkapjuk a gl_LocalInvocationID változóból (ugyanúgy, mint fent a DFT-nél). Ezzel a műveletigény közel log₂N lesz, azaz például egy 1024 hosszú mintasorozatra 10 lépésben (!) kiszámolja a DFT-t. A kódot most még nem írom le (bár nem kell hozzá nagy fantázia), majd a konkrét óceán spektrumok (inverz) transzformációjánál.

(megj.: nem tudom észrevettétek-e, hogy a pillangó műveletből eltűnt u - t, így műveletigényben majdnem ugyanolyan gyors ez is, mint az előző. Höfö rájönni, hogy miért :) Annyit segítek, hogy a twiddle factor-ral kapcsolatos)

A matematikai háttérrel felvértezve térjünk át az óceánfelület modelljeire. Ezt a tudományágat tengertannak vagy fizikai oceanográfiának nevezik. Ebbe beletartozik a hullámok viselkedésének vizsgálata, illetve a fény-óceán interakció is (de utóbbiba nem megyek bele részletesen).

Először ismertetek néhány alapfogalmat, amikre majd építkezni fogok. A folyadékdinamika Euler-féle egyenleteinek egy pontos megoldását (ún. Gerstner hullámok) már kb. 200 éve felfedezték, ami pongyolán a következőt mondja: vegyünk egy x₀ pontot és egy y₀ = 0 magasságot a nyugalomban levő óceánfelületen. A Gerstner modell azt írja le, hogy amikor egy k irányú, A amplitúdójú hullám áthalad a ponton, akkor ezek az értékek hogyan változnak meg. Az észrevétel az, hogy a pont körszerű mozgást végez:

Ahol k a hullámvektor hossza (kapcsolata a hullámhosszal: k = 2π / λ), t az idő, ω pedig a (kör)frekvencia, aminek a kapcsolatát a hullámvektorral az ún. diszperziós reláció írja le.

Ez láthatóan egy viszonylag egyszerű modell, nincs is túl sok köze a valósághoz, viszont több ilyen összekombinálásával elő lehet állítani többé-kevésbé hihető vízfelületet, például ha a célhardver nem elég erős a komplikáltabb módszerekhez.

A diszperziós reláció más lehet attól függően, hogy milyen vízfelületről van szó. Két dolog befolyásolja: a víz mélysége (D), illetve a felületi feszültség (ez most nem érdekes). Az általános képlet


ahol g = 9.81 m/s² a gravitációs konstans. Amennyiben a szóban forgó víz elég mély, akkor a hiperbolikus tangens eltűnik (gondolatban felrajzolva kD = 3-tól már nyugodtan tekinthető 1-nek). A megjelenítés szempontjából is célszerű mély óceánnal foglalkozni, mert akkor nincs szükség fénytörés számolásra (amiről tudjuk, hogy problematikus ilyen jellegű felületeknél).

Dióhéjban a Gerstner modellel az a gond, hogy nem trivi elérni vele, hogy periodikus legyen (megoldható, de mivel eleve nem realisztikus, fölösleges vacakolni). Persze ez inkább 3D-ben jelent problémát, 2D-ben nyugodtan lehet használni, sőt fel is adom höfönek, hogy csináljatok egy egyszerű 2D játékot, amiben ilyen hullámok vannak.

---

Addig is ez kuka, térjünk át a statisztikai alapú modellre. Erről azt kell tudni, hogy az óceánfelület magasságát valószínűségi változóként kezeli, az oceanográfiából származó eredmények pedig kimutatták, hogy ezek a magasságok következetesen a normális (Gauss) eloszlást követik.

Az analitikus modellek építéséhez szükséges adatokat a hullámok magasságának mintavételezésével szerzik meg, pl. bójákkal, radarral vagy egyszerűen csak megfigyeléssel. Ha valakit mélyebben érdekel a mérés menete, illetve az adatokból előálló spektrum, akkor itt egy rövid összefoglaló.

Most viszont analitikus modellt keresünk, ami lehetőleg minél jobban közelíti a mért adatokat. Ez hasonlóan mint előbb, az óceánfelület magasságának szinusz illetve koszinusz függvényekre bontásán alapul, és az alábbi képlet írja le (ami láthatóan egy kétdimenziós inverz DFT):


Ahol (előrevetítve, hogy egy displacement map-be lesz majd kiszámolva)


és N, M a displacement map felbontása, L_x, L_z pedig az óceándarakba térbeli méretei (méterben). A fenti kétdimenziós inverz DFT az x = (nL_x/N, mL_z/M) diszkrét pontokban adja meg a hullámok magasságát (ezt később ki fogom használni).

(megj.: ezeket "lineáris" vagy "gravitációs" hullámoknak nevezik, de sekély vízben megjelennek "nemlineáris" hullámok is)

A következő lépés egy analitikus hullám spektrum (jelölésben P_h(k)) kiválasztása, amiből a DFT együtthatói előállíthatóak. Két ilyen nevezetes spektrumot említek meg, az egyik a Pierson-Moskowitz spektrum, a másik a Phillips spektrum (erről nem találtam irodalmat), de mivel utóbbi az elterjedtebb, azt fogom használni. A képlete


(megj.: HTML-ben nem tudok kalapot írni, itt normalizálást jelent)

ahol L = V² / g a legmagasabb hullám egy adott konstans V szélerősség mellett, w a szintén konstans szélirány, A pedig egy megválasztható konstans, amivel a hullámok magasságát lehet változtatni (általában egy kicsi érték).

Tessendorf szerint a gyorsabb konvergencia érdekében célszerű az L-nél lényegesen kisebb (ε) hullámokat elhagyni, ami a képlet exp(-k²ε²)-tel való megszorzását jelenti (erre innentől módosított Phillips spektrumként hivatkozok).

Az előbbi módosított Phillips spektrumot véletlenszerű értékekkel megszorozva, a Fourier együtthatókra (amplitúdókra) az alábbi képlet használható a t = 0 időpontban

ahol ξ_r és ξ_i egy 0 várható értékű, 1 szórásnégyzetű normális eloszlásból húzott véletlenszerű értékek (ehhez lehet használni az std::normal_distribution osztályt, vagy a Box-Muller módszert). Ennek és a frekvenciáknak (ld. diszperziós reláció) a kiszámolását rátok bízom; csak le kell kódolni a képleteket (az eredményt pedig float textúrákba pakolni).

Ha ezek megvannak, akkor a Fourier együtthatókat a következőképpen lehet kiszámolni valamilyen tetszőleges t időpontban:


Mint később látni fogjuk, a kiinduló állapot szimmetriája miatt nem kell külön kiszámolni a negatív hullámvektorhoz tartozó együtthatókat. Erre a képletre innentől magasságspektrumként fogok hivatkozni.

Ha elvégezzük rá az inverz DFT-t, akkor egy majdnem hihető óceánfelület-darabkát kapunk eredményül. A helyzet viszont az, hogy az óceánfelület nem csak a magasságában (heightfield) fluktuál, hanem az alapsíkjában is, ezáltal fodrozódásokat keltve (choppy field). Megfelelő magasságú hullám esetén ezek a fodrozódások egy hurkot állítanak elő a felületben, ettől ott az óceán habzani fog. Nézzük meg először a különbséget:


A képeken pontosan ugyanabban az időpontban látható egy csak magasságban, illetve egy magasságban és alapsíkban is hullámzó felület. Utóbbi láthatóan "érdesebb", és jobban is emlékeztet az óceánra. Ennek eléréséhez az ún. Lie transzformáció 2D-s megoldását extrapolálja a dolgozat 3D-be. Az alapsíkbeli elmozdulás képletére így a következőt kapjuk:


Ami szintén kétdimenziós inverz DFT, de szerencsére a Fourier együtthatói a magasságspektrumból kiszámolhatóak, ezáltal plusz két spektrumot kapva (gondolom mindenki látja, hogy vastag betűvel írtam, merthogy vektorértékű). A DFT linearitását kihasználva ez a két spektrum összevonható egybe.

Mindegyik spektrumot (a kiindulásit is beleértve) kétdimenziós RG32F textúrákban tárolom, ami az ún. prefetch cache miatt hatékony compute shader implementációt tesz lehetővé. A kód kicsit hosszú, úgyhogy rövidíteni fogom:


Amint látható egy az egyben átírtam kódba a fenti képleteket; ennek fényében gondolom nem kell elmagyaráznom, hogy a kiindulási spektrumot hogyan kell kiszámolni. Annyit fűznék csak hozzá, hogy a hullámvektorokat megfordítottam; nem jöttem rá ugyanis, hogy az eredményül kapott óceán miért ellentétes irányba hullámzik, mint kellene (a legjobb tippem az, hogy az OpenGL fordított koordinátarendszerei miatt). A kód által előállított spektrumok (skálázás után) a következőképpen festenek:


A sorrend ugyanaz, mint ahogy bemutattam; az első tehát a magasságspektrum, a másik kettő pedig az x és z-beli elmozdulás spektrumai (szétszedve). Mindegyiknél megfigyelhető az átlókra való szimmetria (ne tévesszen meg a spektrum alakja!), emiatt optimalizációs célból elég lenne csak az egyik felüket kiszámolni, de az megnehezítené a textúra címzést.

(megj.: próbáljátok ki a kódmellékletben szétszedés nélkül; meg fogtok lepődni :) a matematikai magyarázata höfö)

Végre elérkeztünk a cikk tetőpontjához, ami nem más, mint az imént kiszámolt spektrumok (kétdimenziós) inverz gyors Fourier transzformációja. Mielőtt ennek nekiállnánk, fejtsük ki a magasságmező inverz DFT-jét, hogy jobban lássuk mit is kell majd csinálni. Ehhez először is felteszem most már, hogy M = N, illetve L_x = L_z = N (itt használom ki az x korábbi felírását). Az előzőleg bemutatott FFT implementáció értelmében a hullámvektort is át kell előbb definiálni:

Mégegyszer azért, mert egyelőre nem térben akarom megkapni a hullámokat, hanem egy N×N-es textúrában (a spektrumokban van a valódi információ!). Ekkor kifejtve a magasságmező képletét (az olvashatóság miatt a fizikában megszokott jelölést használva) az alábbit kapjuk:


A kiemelések a már említett e^iπ = -1 miatt tehetőek meg (ha nem hiszitek, ellenőrizzétek az Euler formulával), az átrendezések pedig az exponenciális függvény azonosságai miatt. A képletből egy különösen fontos dolog látszik: a kétdimenziós DFT szeparálható (ugyanúgy, mint pl. a Gauss szűrő). Elég tehát egydimenziós esetre megírni az inverz FFT algoritmust, majd lefuttatni először a spektrum soraira, utána pedig az oszlopaira. Lássuk akkor most már a teljesértékű compute shader-t:


Ravasz módon úgy írtam meg, hogy ne kelljen az indexeléssel bajlódni (viszont értelemszerűen kell egy textúra a köztes adatnak). A képletből kiemelt előjelet az érthetőség kedvéért ide raktam, de természetesen megspórolható az elágazás, ha nem itt végzem el, hanem a három hullámtextúra egyesítésekor (a kódmellékletben természetesen úgy szerepel).

A workgroup shared memória méretkorláta miatt ez az algoritmus legfeljebb 1024×1024-es textúrákra működik, de ez annyira nem nagy baj, mert efölött már látható precíziós hibák jönnek elő. Ha esetleg nagyobb méretre van szükség, akkor globális memóriába pakolt bufferekkel az is megoldható, viszont értelemszerűen lassabb lesz.

Na de lássuk az eredményt (mindhárom spektrumra):


A képeken nem látszik, de higgyétek el, hogy mindegyik a neki megfelelő irányba mozog, illetve az amplitúdóknak megfelelően változik. Ezt a hármat összepakolva egy RGBA32F textúrába elő is áll a displacement map, amivel már lehet játszadozni. Ha eddig nem lett volna világos, a DFT miatt természetesen mindegyik periodikus.

A hatékony de mégis látványos megvilágításhoz szükség van egy normal map-re, illetve az említett "habzáshoz" egy folding map kiszámolására (utóbbit az nVidia nem használja semmire, de majd megmutatom egy lehetséges alkalmazását).

A normálvektorok kiszámítására két lehetőség is van: az egyik, hogy kiszámoljuk a h(x, t) magasságmező gradiensét (∇), ez azonban egy újabb FFT-t jelentene. A másik lehetőség a véges differencia módszer, ami annyit jelent, hogy a displacement map minden texelére egy 2×2-es környezetben vesszük az eltéréseket (szóval majdnem ugyanaz, mint amit a dFdx és dFdy függvény csinálna).

(megj.: L alatt innentől L_x = L_z-t értem)

Két okból nem használtam az említett GLSL függvényeket: egyrészt mert csak a fragment shader-ben érhetőek el, másrészt a kép széleinél helytelen eredményt adnak. Itt említem meg, hogy a % N és a & (N - 1) között az a lényeges különbség, hogy előbbi megtartja az előjelet, míg a másik átfordul (és egy periodikus textúra esetén nyilván ez a cél).

A fodrozódás kiszámolásához hasonló dolgot kell tenni az ún. Jacobi mátrix formájában, ami szintén parciális deriváltakat tartalmaz, viszont vektorértékű függvényekre (és a fenti choppy field ilyen). A mátrix valamilyen x pontból az x + λD(x, t) pontba való transzformációt írja le, ahol λ egy megválasztható konstans, amivel növelni/csökkenteni lehet a kilengéseket (én fixen λ = 1.3-at használok):


(megj.: a Cauchy-Riemann egyenletek miatt ez megegyezik az f = D_x + iD_y függvény komplex deriváltjával)

A mátrix determinánsa mondja meg a transzformáció egyediségét: amíg nem metszi önmagát (a hullám nem borult át magán), addig egyedi, invertálható, a determináns pedig pozitív értékű. Ennek kiszámolása szintén véges differenciákkal történik:


Gondolom a 2×2-es mátrix determinánsának kiszámolására csak emlékeztek... Zárójelben, amikor a mátrix négyzetes, akkor az irodalmak egyszerűen csak Jacobian-nak nevezik (beleértve a determinánst is).


A gradienseket a szokásos normalize/mad függvényekkel alakítottam a megszokott formára, de ne felejtsük el, hogy az óceán (kvázi) sík felület, ami miatt nem kell tangent frame számolásokkal küzdeni. A folding map használatáról később írok.

Most egy olyan rész jön, amit csak nagyvonalakban mutatok be, ugyanis ezt a kódot speciel egyszerűbb kimásolni, mint oldalakon keresztül magyarázni. Apró módosításokat végeztem csak rajta, például javítottam az indexelést, eltüntettem a degenerált háromszögeket, stb. Az egész egyébként az ún. geomipmapping elméleten alapul.

Kétféle megoldás van, az egyik egy screen space-ben adott háromszöghálót visszatranszformál world space-be, elvégzi az eltolásokat, majd visszavetíti screen space-be. Az előnye, hogy kevés memóriával gyakorlatilag tetszőlegesen messzire ki tudja húzni az óceánfelületet. Hátránya, hogy már a kamera közvetlen közelében is azonnal lecsökken a részletessége. Játékok általában ezt alkalmazzák (pl. mindhárom Crysis), mivel eleve kis felbontású (64×64) textúrákkal dolgoznak (eml.: akkoriban még nem volt compute shader!).

A másik megoldás ezek szerint world space-ben dolgozik, viszont sokkal szofisztikáltabb algoritmusok kellenek ahhoz, hogy elég nagy területet le tudjon fedni. Bár mindkét megoldáshoz találtam példakódot, inkább ezt választottam, mert így én tudom befolyásolni, hogy mikortól kezdjen el csökkenni a részletesség.

A módszer a következő: egy adott óceándarabkát (patch) két részre oszt; egy belső triangle strip hálóra, illetve egy külső, más LOD (level of detail) szintekkel való kapcsolódási pontokat tartalmazó triangle list-re. Ki lehet számolni, hogy ha például három LOD szinttel való összekapcsolódást szeretnénk, akkor egy adott patch-nek összesen 81 változata lehet. Mondjuk 6 LOD szint esetén, akkor ez 486 darab teljesen különböző patch-et jelent. Ezeknek a generálását nem részletezem.

(megj.: a vertexek közötti távolság lehetőleg 2 cm-nél nagyobb legyen; ez alatt ugyanis már előjönnek az említett "nemlineáris" folyamatok, illetve a felületi feszültség)

Feltételezve, hogy a teljes óceán mondjuk 25 km², egy patch pedig 400 m², könnyen rá lehet jönni, hogy valamiféle gyorsítóstruktúra nélkül pillanatok alatt agyonvágja a GPU-t. Ez a gyorsítóstruktúra egy quadtree, amit hála az égnek nem fordítottak le magyarra (javaslom a kvartáris fa elnevezést); a lényeg, hogy egy adott csúcsnak négy darab alcsúcsa van. A gyökércsúcs természetesen a teljes óceánfelület, amiből rekurzív négy részre osztással áll elő a teljes fa. Azaz, ha el akarjuk dönteni, hogy egy adott patch látszik-e, akkor az log₄ 65536 = 8 lépésben megtehető.

Persze elég nagy ostobaság egy ilyen bazinagy fát állandóan benntartani a memóriában, ezért kihasználva, hogy a patch-ek maguk is egy hálót alkotnak, a view frustum-ba eső patch-eket minden rajzoláskor meg lehet határozni anélkül, hogy magát a fát eltárolnánk. Ez tipikusan 130 darab patch-et jelent, ami sokkal barátságosabb szám. A szükséges LOD szintet meg lehet határozni távolságból vagy lefedettségből (coverage; hány pixelt fed le a patch). Kívülről nézve az eredmény valami ilyesmi:

Ezen a képen persze látszik, hogy hol történnek a LOD szint váltások, de a generált kapcsolódási pontok miatt a konkrét óceánon egyáltalán nem vehetőek észre. Lehetne egyébként használni tessellation shader-t is, de sajnos a tesszellációs szint le van korlátozva 64-re, ami nekem nem volt elég (de az ARM/Mali mutat példát).

Korábban lerögzítettem, hogy csak mély óceánnal foglalkozok, azzal is csak felülről (ha lemegyek alá, akkor megint egy BSDF modellt kellene alkalmazni). Ha valakit részletesen érdekel a sekély víz, illetve óceán alatti effektek, akkor javaslom a Crysis ingyenesen letölthető demóját, amiben az összes ezzel kapcsolatos HLSL shader megtalálható.

Kezdjük azzal, hogy a mély óceán felülete szinte tökéletes tükröző (hasonlóan, mint a fémek), azaz a megjelenésének egy jelentős része tükröződésekből áll. A Fresnel egyenletekben tehát figyelmen kívül hagyható a fénytörés; elég valamilyen konstans színt használni helyette. Egy kevésbé szimpatikus tulajdonsága, hogy anizotróp, ami olyan felületet jelent, amit a tengelye körül forgatva változik a kinézete (pl. bársony, haj, csiszolt fémfelületek, stb.).

A víz vagy vizes felületek esetén az anizotropitás abban mutatkozik meg, hogy a tükrözött fények "megnnyúlnak". Ezt leginkább este lehet megfigyelni például a Duna felszínén, vagy esős időben autózva. Két konyhanyelvi leírást ajánlok, amik való életbeli képeken keresztül mutatják be ezt az effektet: egyik, másik.

A fizikai magyarázata ennek a jelenségnek az, hogy a fényforrás egy konkrét pontja több mikrofelületen is tükröződik. Ebből rögtön adódik, hogy a Torrance-Sparrow mikrofelület modellnek megfelelő BRDF-et kellene találni; természetesen olyat, ami az anizotropitást is figyelembe veszi. Előrevetítem, hogy ezek a modellek bonyolultak. Van azonban egy kerülőút, ugyanis a Blinn-Phong modell is bizonyos szintig kezeli ezt a dolgot, de sajnos nem fizikai alapú (a cikk végén levő képeken látszik majd a különbség).

Több ilyen BRDF modell is van, az egyik a Ross-féle (ez nagyon realisztikus), a másik a Ward modell. Mivel utóbbi a könnyebbik, azt implementáltam le, viszont a cikk utolsó simításainak fázisában, szóval lehet, hogy néhány dolgot elrontottam benne. A képlete (megint fizikai jelöléssel mert nem fér ki):

Ahol ρ_s a BRDF lufi "hosszát" jelenti, α_x és α_y az x illetve y iránybeli "vastagságát", h pedig a normalizálatlan félvektor. A teljes shader kód az alábbi módon fest, de a Ward modellhez saját kútfőből számoltam ki az értékeket; nem kerestem referenciát:


Az értékek változtatásával lehet erősíteni, szűkíteni, nyújtani a tükrözött fényt. Annyit elmondanék, hogy a demóban nem HDR-ben rajzolok (nem találtam megfelelő égboltot), a Ward modell 1-hez való integrálását pedig nem ellenőriztem, amitől lehet, hogy egy kicsit "erős" (ezen a paraméterek buherálásával lehet segíteni). A szélességét szintén empirikusan határoztam meg, hogy az égbolthoz nagyjából igazodjon (de a valóságban ennél jóval keskenyebb).


Nyilván egy teljesértékű óceánhoz ez nem elég. A Ross modellhez itt található egy mindent részletező cikk, amiket én nem csináltam meg (például sík tükröződéseket). Ha valaki késztetést érez rá, akkor javaslom, hogy screen space reflection módszerrel próbálkozzon.

Utolsó erőfeszítésként az nVidia javaslatai alapján tüntessük el a periodicitásból adódó (és igen zavaró) ismétlődéseket. Az ötlet annyi, hogy egy bizonyos távolságtól kezdve cseréljük le az FFT alapú óceánt egy Perlin zaj alapú hullámra, ami érdekes módon még ismételve sem produkál látható artifact-okat (így például használják domborzat, felhő, stb. megjelenítéséhez is). A képlete leegyszerűsítve (saját formalizáció):

Ahol P az oktávok száma, A_j és f_j a j-edik oktáv amplitúdója illetve frekvenciája, p(x) a vetített textúra koordináta, noise pedig az alap Perlin textúra. A hullám, illetve a normálvektorai közvetlenül generálhatóak shader-ben:




A többi pedig ugyanaz, mint eddig. Összehasonlító kép:


A módszer akkor is működik, ha elég magasról nézek közvetlenül az óceánfelületre, ekkor viszont a kép közepén lesz részletes és a szélei felé vált át Perlin hullámokba.

Az ehhez hasonló mély matematikai cikkeknél általánosan jellemző, hogy az implementáció gyorsan megvan, az elméleti háttér szöveggé formálása viszont rengeteg időt elvisz. Adódik ez például abból, hogy egy adott matematikai fogalom felvezetése után bennem is felmerülnek kérdések, hogy "na de miért van ez így?". A válasz felkutatása során pedig még mélyebbre ásom magam a már egyébként is nehéz elméletekben.

A cikk eredeti célja az volt, hogy megtanuljam a Fourier transzformációt (mert mint említettem, kimaradt a tanulmányaimból). Mivel tudtam azonban, hogy az óceán megjelenítése a Waterworld (1995) illetve a Titanic (1997) óta ezen az elméleten alapul, egyértelművé vált, hogy így fogom bemutatni. A végeredmény viszont a becslésemmel ellentétben erősen a matematikai oldal felé húz, és kevésbé fókuszál az óceán fizikai alapú megjelenítésére (bár természetesen a megfelelő hivatkozásokat megadtam).

Hová lehetne akkor ebben a témakörben tovább kutatni? Említettem például a folyadékdinamikát (azon belül a Navier-Stokes egyenleteket); ez inkább fizika, annyira nem áll hozzám közel. Felmerült az óceán (fény)fizikailag teljesértékű megjelenítése (példaként a Ross-féle BRDF modellt említettem meg); ez külön cikket érdemelne.

Ezek mellett a bevezetésben hivakoztam a 20. század legjelentősebb algoritmusaira, amik közül a legelső természetesen a pónilósimogató AI (csak viccelek; az első a Monte Carlo Metropolis Light Transport), amin én is meglepődtem. Akármelyik irányba is megyek majd tovább, az a léc egy újabb megemelését jelenti majd. Az alábbi képeken kétféle alapszínnel mutatom be az eredményt:

(megj.: bal oldal: Ward modell, jobb oldal: Blinn-Phong modell; előbbi a valóságnak megfelelve elér teljesen az óceán végéig)

Addig is a kódmelléklet elérhető a szokott helyen, a teljes mértékig felturbózott óceánról készült videó pedig itt tekinthető meg (állítsátok át 1080p-re, ha nem abban indul). A diagramok gnuplot kódját kérésre odaadom.


Höfö:

• Csinálj egy 2D játékot Gerstner hullámokkal!
• Cseréld le a Phillips spektrumot a Pierson-Moskowitz-re!
• Implementáld le a screen space-beli megoldást!

Irodalomjegyzék

http://evasion.imag.fr/.../Jonathan/articlesCG/waterslides2001.pdf - Simulating Ocean Surfaces (Jerry Tessendorf, 2001)
https://people.cs.clemson.edu/~jtessen/papers_files/coursenotes2004.pdf - Simulating Ocean Water (Jerry Tessendorf, 2004)
https://pdfs.semanticscholar.org/330e/... - Notes on the Ward BRDF (Bruce Walter, 2005)
http://www-evasion.imag.fr/~Fabrice.Neyret/.../NV_OceanCS_Slides.pdf - Ocean Surface Simulation (nVidia)

https://hal.inria.fr/inria-00443630/file/article-1.pdf - Real-time Realistic Ocean Lighting [...] to BRDF (EUROGRAPHICS, 2010)
https://dspace5.zcu.cz/bitstream/11025/11943/1/Tian.pdf - Ocean wave simulation by the mix of FFT and Perlin (WSCG, 2014)
https://.../Realtime_GPU_ocean_rendering.pdf?dl=0 - Real-time GPU-driven Ocean Rendering on Mobile (ARM/Mali, 2015)
http://developer.download.nvidia.com/... - Ocean simulation and rendering in War Thunder (nVidia, 2015)
https://tubdok.tub.tuhh.de/... - Realtime GPGPU FFT Ocean Water Simulation (Hamburg University of Technology, 2017)

http://mogi.bme.hu/TAMOP/mereselmelet/ - Méréselmélet (BME, 2014)