Spotlight #62 – WebGL 3D — GLSL Raymarching, Pola Odległości ze Znakiem i Symulacja GPGPU

I. Raymarching i Pola Odległości ze Znakiem — Podstawowa Technika

Zanim przyjrzymy się poszczególnym symulacjom, warto zrozumieć wspólny fundament, na którym wszystkie się opierają. Tradycyjne renderowanie 3D działa poprzez przekształcanie trójkątów w potoku rasteryzacji: każda powierzchnia musi zostać opisana jako siatka wielokątów. Sprawdza się to znakomicie dla modeli inżynierskich, ale słabo pasuje do obiektów zdefiniowanych niejawnie — fraktali, powierzchni płynów, kształtów algebraicznych — gdzie żadna siatka trójkątów nie jest praktyczna.

Pole odległości ze znakiem (SDF) to funkcja d(p), która zwraca odległość ze znakiem od dowolnego punktu p w przestrzeni 3D do najbliższej powierzchni: dodatnią na zewnątrz, ujemną wewnątrz, zero na powierzchni. Dla kuli wyśrodkowanej w początku układu współrzędnych o promieniu r:

d(p)  =  length(p)  -  r

Pola SDF można łączyć prostą arytmetyką: sumą dwóch kształtów jest min(d1, d2); przecięciem jest max(d1, d2); odejmowanie (wydrążenie jednego kształtu z drugiego) to max(d1, -d2). Ta algebra kształtów, zwana konstruktywną geometrią bryłową (CSG), pozwala budować dowolnie złożone obiekty z prymitywów bez konieczności konstruowania siatki.

Pętla Raymarchingu

Mając SDF, raymarching renderuje go, wystrzeliwując promień z kamery przez każdy piksel ekranu i maszerując wzdłuż promienia krokami o wielkości wyznaczonej przez wartość SDF:

float t = 0.0;
for (int i = 0; i < MAX_STEPS; i++) {
    vec3 p = ray_origin + t * ray_dir;
    float d = sdf(p);
    if (d < EPSILON) { /* trafienie w powierzchnię */ break; }
    if (t > MAX_DIST) { /* pudło, kolor nieba */ break; }
    t += d;  // bezpieczny krok: nie przekroczy powierzchni
}

Kluczowa intuicja polega na tym, że wartość SDF jest bezpiecznym rozmiarem kroku: promień nie może przekroczyć powierzchni, ponieważ d(p) jest dokładnie minimalną odległością do dowolnej powierzchni z punktu p. Ten algorytm sphere tracing (Hart, 1996) jest zarówno bezpieczny, jak i wydajny: w otwartej przestrzeni wykonywane są duże kroki; w pobliżu powierzchni kroki kurczą się do precyzji subpikselowej. Typowa złożona scena zbiega się w 50–100 iteracjach na piksel, każda iteracja wykonuje jedną ewaluację SDF — obciążenie idealnie pasujące do masowo równoległego wykonania na GPU, z jednym wątkiem na piksel.

Normalne Powierzchni i Oświetlenie

Gdy trafienie zostanie znalezione, normalna powierzchni jest obliczana metodą różnic skończonych — próbkując SDF w sześciu pobliskich punktach i biorąc gradient:

vec3 normal(vec3 p) {
    float e = 0.001;
    return normalize(vec3(
        sdf(p + vec3(e,0,0)) - sdf(p - vec3(e,0,0)),
        sdf(p + vec3(0,e,0)) - sdf(p - vec3(0,e,0)),
        sdf(p + vec3(0,0,e)) - sdf(p - vec3(0,0,e))
    ));
}

Wymaga to tylko sześciu dodatkowych wywołań SDF na piksel i działa dla dowolnej powierzchni niejawnej bez jawnego obliczania normalnej siatki. Standardowe oświetlenie Blinna-Phonga lub fizycznie oparte wykorzystuje następnie normalną, kierunek widzenia i kierunek światła w zwykły sposób. Miękkie cienie i okluzja otoczenia są obliczane przez dodatkowe krótkie marsze promienia w stronę źródła światła: jeśli minimalna wartość SDF promienia cienia jest bardzo mała, punkt jest częściowo zasłonięty, co daje gładkie półcienie niemożliwe przy twardych mapach cieni.

II. Mandelbulb — Fraktal 3D z Estymacją Odległości

Zbiór Mandelbrota żyje na płaszczyźnie zespolonej: dla każdego punktu c iteruj z → z² + c i sprawdzaj, czy orbita ucieka. Rozszerzenie tego na trzy wymiary wymaga pojęcia „mnożenia" dla wektorów 3D. Mandelbulb (White i Nylander, 2009) używa potęgi triplex we współrzędnych sferycznych (r, θ, φ):

z^n:  r_new = r^n
      θ_new = n · θ
      φ_new = n · φ

Nie jest to algebra w ścisłym sensie (brakuje jej łączności dla n ≠ 2), ale daje wizualnie niezwykłe rezultaty dla n = 8: bulwiasta struktura pokryta zawiłymi sferycznymi wybrzuszeniami, z których każde, przy bliższym oglądzie, pokryte jest mniejszymi kopiami tych samych wybrzuszeń — samopodobny fraktal na wszystkich skalach.

Estymator Odległości

Mandelbulb nie ma zamkniętej formy SDF, ale ma estymator odległości (DE) wyprowadzony z pochodnej mapy iteracji. Śledź zarówno bieżącą iterację z, jak i jej pochodną dz/dc podczas pętli iteracji:

// Każda iteracja:
dz = n * pow(r, n-1) * dz + 1.0;
z  = triplex_pow(z, n) + c;
if (length(z) > BAILOUT) break;

Po ucieczce iteracji estymata odległości wynosi:

de  =  0.5 · log(|z|) · |z| / |dz/dc|

Ta estymata jest dokładna z dokładnością do małego stałego czynnika w pobliżu powierzchni, co wystarcza do raymarchingu: pętla sphere tracing może bezpiecznie użyć jej jako rozmiaru kroku. Kosztem jest to, że każdy krok raymarchingu wymaga uruchomienia pełnej iteracji fraktala, by wyznaczyć DE, co czyni Mandelbulba jedną z najbardziej obciążających GPU symulacji na stronie — a jednak wciąż interaktywną przy 60 fps na nowoczesnym sprzęcie.

III. Gąbka Mengera — IFS przez Składanie SDF

Gąbka Mengera powstaje przez wielokrotne usuwanie centralnej kostki i centrów każdej ściany z sześcianu — iterowany system funkcji (IFS), który daje fraktal o wymiarze Hausdorffa log(20)/log(3) ≈ 2,727. Zbudowanie siatki nawet dla 4 iteracji wymaga milionów wielokątów. Podejście SDF, wykorzystujące składanie dziedziny, ewaluuje dowolną głębokość iteracji funkcją o stałym koszcie.

Kluczowym trikiem jest składanie modularne: zamiast rekurencyjnie dzielić dziedzinę, składamy przestrzeń tak, by pojedyncza ewaluacja SDF podstawowej kostki reprezentowała cały fraktal. Przy każdej iteracji współrzędne są składane do fundamentalnej dziedziny grupy symetrii, skalowane o czynnik 3, a krzyżowa wnęka jest odejmowana:

float menger_sdf(vec3 p, int iterations) {
    float d = box_sdf(p, vec3(1.0));
    float scale = 1.0;
    for (int i = 0; i < iterations; i++) {
        // Złóż do komórki jednostkowej, odbij wzdłuż osi
        vec3 q = mod(p * scale, 2.0) - 1.0;
        q = abs(q);
        // Wnęka krzyżowa: usuń centralną rurę w każdej parze osi
        float cross = min(max(q.x, q.y),
                     min(max(q.y, q.z),
                         max(q.z, q.x)));
        float cavity = cross - 1.0/3.0;
        d = max(d, -cavity / scale);
        scale *= 3.0;
    }
    return d;
}

Każda iteracja pętli odpowiada jednemu poziomowi hierarchii fraktala. Dzielenie przez scale koryguje SDF, uwzględniając kompresję przestrzenną wprowadzoną przez złożenie. Piękno tego podejścia polega na tym, że iteracje od 1 do 6+ są dostępne kosztem kilku dodatkowych iteracji pętli, bez zmiany geometrii — tylko suwak kontrolujący liczbę iteracji.

IV. Wormhole 3D — Geometria Morrisa-Thorne'a w GLSL

Tunel czasoprzestrzenny Morrisa-Thorne'a jest rozwiązaniem równań pola Einsteina opisującym przejezdny tunel między dwoma regionami czasoprzestrzeni. Metryka we współrzędnych sferycznych to:

ds²  =  -e^(2Φ) dt²  +  dl²  +  (b₀² + l²)(dθ² + sin²θ dφ²)

gdzie l to odległość własna wzdłuż osi tunelu, b₀ to promień gardła (minimalny promień obwodu / 2π), a Φ(l) to funkcja przesunięcia ku czerwieni. Powierzchnia osadzenia — dwuwymiarowy przekrój oddający geometrię przestrzenną — to powierzchnia obrotu z(r) = ±b₀ · ln(r/b₀ + √((r/b₀)²-1)), która rozszerza się na zewnątrz od gardła jak dwie połączone trąbki.

Symulacja GLSL nie rozwiązuje pełnych ogólnorelatywistycznych równań promieni (co wymagałoby całkowania sprzężonych równań różniczkowych na piksel na klatkę). Zamiast tego wykorzystuje przybliżenie geometryczne: promienie są zaginane w miarę zbliżania się do gardła przez shader wypaczający przestrzeń, który płynnie mapuje kierunki wejściowe na wyjściowe zgodnie z wewnętrzną krzywizną osadzenia. Gwiazdy i mgławice zdają się owijać wokół gardła w wyniku soczewkowania grawitacyjnego. To ta sama technika, użyta w filmie Interstellar (2014), dla którego Kip Thorne i zespół efektów wizualnych Double Negative opracowali fizycznie dokładny renderer GLSL tunelu czasoprzestrzennego.

V. Supershape — Superformuła Gielisa na GPU

Johan Gielis wprowadził superformułę w 2003 roku jako uogólnienie okręgu, które potrafi opisać niezwykłą różnorodność naturalnych i matematycznych kształtów zaledwie sześcioma parametrami (m, n1, n2, n3, a, b):

r(θ)  =  [ |cos(mθ/4)/a|^n2  +  |sin(mθ/4)/b|^n3 ]^(-1/n1)

We współrzędnych biegunowych r(θ) daje promień pod kątem θ. Dla supershape 3D dwie niezależne ewaluacje superformuły są stosowane jako promienie szerokości i długości geograficznej, dając powierzchnię przypominającą kulę zdeformowaną przez iloczyn dwóch funkcji promieniowych. GPU mapuje każdy piksel bezpośrednio na współrzędną (θ, φ), ewaluuje dwie formuły, oblicza pozycję 3D na powierzchni i renderuje przez standardowy model oświetlenia — bez siatki, bez shadera geometrii, tylko czysta powierzchnia parametryczna ewaluowana na piksel.

VI. Cząstki GPGPU — 250 000 Cząstek przez Render-to-Texture

Centralnym wyzwaniem symulacji cząstek na GPU jest ograniczenie WebGL: shader nie może jednocześnie czytać z i pisać do tej samej tekstury. Technika ping-pong rozwiązuje to, utrzymując dwie pary celów renderowania (A i B). W nieparzystych klatkach przejście obliczeniowe czyta z A i pisze do B; w parzystych czyta z B i pisze do A. CPU nigdy nie dotyka danych cząstek po inicjalizacji.

// Pseudokod dla jednej klatki:
gl.bindFramebuffer(writeFBO);     // powiąż cel B
gl.useProgram(computeShader);
gl.bindTexture(readTexture);      // powiąż cel A
gl.drawArrays(TRIANGLES, 0, 6);  // pełnoekranowy quad: uruchamia fizykę na GPU
swap(A, B);                       // następna klatka czyta z B

// Przejście renderowania:
gl.bindFramebuffer(null);         // bufor tylny
gl.useProgram(renderShader);
gl.bindTexture(positionTexture);  // bieżące pozycje (właśnie zapisane)
gl.drawArrays(POINTS, 0, N);     // każdy wierzchołek punktu próbkuje własny teksel

Pozycja każdej cząstki jest przechowywana jako teksel RGBA32F (cztery 32-bitowe liczby zmiennoprzecinkowe: pozycja xyz, w dla wieku lub masy). Tekstura 512×512 przechowuje 262 144 cząstek. Shader obliczeniowy to pełnoekranowy quad, którego shader fragmentów wyszukuje pozycję cząstki w gl_FragCoord, stosuje aktualizację fizyki (grawitacja, szum curl lub atraktor) i zapisuje nową pozycję jako kolor wyjściowy. Shader renderujący to wywołanie rysowania punktów: każdy indeks wierzchołka punktu mapuje się na współrzędną teksela, odczytuje pozycję przez pobranie tekstury w shaderze wierzchołków i wyprowadza sprite punktowy.

Szum Curl i Przepływ Bezdywergencyjny

Najbardziej wizualnie przekonujące zachowanie cząstek wykorzystuje szum curl dla pola prędkości. Mając dowolny 3D potencjał wektorowy A(p), curl F = ∇ × A jest matematycznie gwarantowany jako bezdywergencyjny:

∇ · (∇ × A)  =  0   (tożsamość rachunku wektorowego)

Pole bezdywergencyjne nie ma źródeł ani ujść — cząstki ani nie zbijają się w gęste kłębki, ani nie rozlatują się. Zamiast tego tworzą nieściśliwe, wirujące strumienie, które wyglądają uderzająco jak dym, chmury i prądy oceaniczne. Potencjał wektorowy A przyjmowany jest jako 3D pole szumu simplex, a curl obliczany jest analitycznie (gradient szumu jest różniczkowalny). Adwekcja 250 000 cząstek przez szum curl przy 60 fps nie wymaga zaangażowania CPU po skompilowaniu shadera — GPU obsługuje cały potok symulacji i renderowania.