Spotlight #62 – WebGL 3D — GLSL Raymarching, поля знакових відстаней та GPGPU-симуляція

I. Raymarching і поля знакових відстаней — основна техніка

Перш ніж розглядати окремі симуляції, варто зрозуміти спільний фундамент, на якому вони всі побудовані. Традиційний 3D-рендеринг працює, перетворюючи трикутники через конвеєр растеризації: кожну поверхню потрібно описати як сітку полігонів. Це чудово працює для інженерних моделей, але погано підходить для об'єктів, заданих неявно — фракталів, поверхонь рідин, алгебричних форм — де жодна трикутна сітка непрактична.

Поле знакових відстаней (SDF) — це функція d(p), що повертає знакову відстань від будь-якої точки p у 3D-просторі до найближчої поверхні: додатну зовні, від'ємну всередині, нуль на поверхні. Для сфери з центром у початку координат і радіусом r:

d(p)  =  length(p)  -  r

SDF можна комбінувати простою арифметикою: об'єднання двох форм — це min(d1, d2); перетин — max(d1, d2); віднімання (вирізання однієї форми з іншої) — max(d1, -d2). Ця алгебра форм, яку називають конструктивною твердотільною геометрією (CSG), дозволяє будувати як завгодно складні об'єкти з примітивів, ніколи не конструюючи сітку.

Цикл raymarching

Маючи SDF, raymarching рендерить його, випускаючи промінь від камери крізь кожен піксель екрана й просуваючись уздовж променя кроками, розмір яких визначається значенням 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) { /* влучили в поверхню */ break; }
    if (t > MAX_DIST) { /* промах, колір неба */ break; }
    t += d;  // безпечний крок: не перескочить поверхню
}

Ключове спостереження в тому, що значення SDF є безпечним розміром кроку: промінь не може перескочити поверхню, оскільки d(p) — це саме мінімальна відстань до будь-якої поверхні з точки p. Цей алгоритм сферичного трасування (Харт, 1996) одночасно безпечний і ефективний: у відкритому просторі кроки великі; поблизу поверхні кроки зменшуються до субпіксельної точності. Типова складна сцена сходиться за 50–100 ітерацій на піксель, кожна ітерація виконує одне обчислення SDF — навантаження, що ідеально лягає на масово-паралельне виконання GPU з одним потоком на піксель.

Нормалі поверхні та освітлення

Щойно знайдено влучення, нормаль поверхні обчислюється скінченними різницями — вибіркою SDF у шести сусідніх точках і взяттям градієнта:

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))
    ));
}

Це вимагає лише шести додаткових викликів SDF на піксель і працює для будь-якої неявної поверхні без явного обчислення нормалі сітки. Стандартне освітлення за Блінном-Фонгом або фізично коректне освітлення потім використовує нормаль, напрямок огляду й напрямок світла звичним чином. М'які тіні й оточувальна оклюзія обчислюються додатковими короткими трасуваннями променів у бік джерела світла: якщо мінімальне значення SDF тіньового променя дуже мале, точка частково затінена, що дає плавні напівтіні, неможливі з жорсткими картами тіней.

II. Mandelbulb — 3D-фрактал з оцінкою відстані

Множина Мандельброта живе на комплексній площині: для кожної точки c ітеруємо z → z² + c і перевіряємо, чи орбіта втікає в нескінченність. Розширення цього до трьох вимірів вимагає поняття «множення» для 3D-векторів. Mandelbulb (Уайт і Ніландер, 2009) використовує триплексний степінь у сферичних координатах (r, θ, φ):

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

Це не алгебра в строгому сенсі (їй бракує асоціативності при n ≠ 2), але вона дає візуально надзвичайні результати для n = 8: бульбашкоподібну структуру, вкриту складними сферичними опуклостями, кожна з яких, при ближчому розгляді, вкрита меншими копіями тих самих опуклостей — самоподібний фрактал на всіх масштабах.

Оцінювач відстані

Mandelbulb не має SDF у замкненій формі, але має оцінювач відстані (DE), виведений із похідної ітераційного відображення. Відстежуйте водночас поточну ітерацію z і її похідну dz/dc протягом циклу ітерацій:

// Кожна ітерація:
dz = n * pow(r, n-1) * dz + 1.0;
z  = triplex_pow(z, n) + c;
if (length(z) > BAILOUT) break;

Після того як ітерація втекла в нескінченність, оцінка відстані дорівнює:

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

Ця оцінка точна з точністю до невеликого постійного множника поблизу поверхні, чого достатньо для raymarching: цикл сферичного трасування може безпечно використовувати її як розмір кроку. Ціна в тому, що кожен крок raymarching вимагає запуску повної фрактальної ітерації для визначення DE, що робить Mandelbulb однією з найінтенсивніших для GPU симуляцій на сайті — проте вона все ж інтерактивна на 60 кадрах/с на сучасному залізі.

III. Губка Менгера — IFS через згортання SDF

Губка Менгера будується повторним видаленням центрального кубика й центрів кожної грані з куба — система ітеративних функцій (IFS), що дає фрактал із розмірністю Хаусдорфа log(20)/log(3) ≈ 2.727. Побудова сітки навіть для 4 ітерацій вимагає мільйонів полігонів. Підхід через SDF, використовуючи згортання простору, обчислює будь-яку глибину ітерації функцією з постійною вартістю.

Ключовий трюк — модулярне згортання: замість рекурсивного поділу простору, простір складається так, щоб одне обчислення SDF базового куба представляло весь фрактал. На кожній ітерації координати згортаються у фундаментальну область групи симетрії, масштабуються втричі, і віднімається хрестоподібна порожнина:

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++) {
        // Згорнути в одиничну комірку, відобразити крізь осі
        vec3 q = mod(p * scale, 2.0) - 1.0;
        q = abs(q);
        // Хрестоподібна порожнина: видаляємо центральну трубку в кожній парі осей
        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;
}

Кожна ітерація циклу відповідає одному рівню фрактальної ієрархії. Ділення на scale коригує SDF, враховуючи просторове стиснення, внесене згортанням. Краса цього підходу в тому, що ітерації від 1 до 6+ доступні ціною лише кількох додаткових ітерацій циклу без жодної зміни геометрії — просто повзунок, що керує кількістю ітерацій.

IV. Червоточина 3D — геометрія Морріса-Торна в GLSL

Червоточина Морріса-Торна — розв'язок рівнянь поля Ейнштейна, що описує прохідний тунель між двома областями простору-часу. Метрика у сферичних координатах:

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

де l — власна відстань уздовж осі тунелю, b₀ — радіус горловини (мінімальний радіус кола / 2π), а Φ(l) — функція червоного зсуву. Поверхня вкладення — 2D-зріз, що відображає просторову геометрію — це поверхня обертання z(r) = ±b₀ · ln(r/b₀ + √((r/b₀)²-1)), яка розширюється назовні від горловини, наче два з'єднані рупори.

GLSL-симуляція не розв'язує повні загальнорелятивістські рівняння променів (що вимагало б інтегрування зв'язаних диференціальних рівнянь на кожен піксель у кожному кадрі). Натомість вона використовує геометричне наближення: промені викривляються, наближаючись до горловини, шейдером просторового викривлення, який плавно відображає вхідні напрямки у вихідні відповідно до внутрішньої кривини вкладення. Зорі й туманності, здається, огортають горловину через гравітаційне лінзування. Це та сама техніка, що використана у фільмі «Інтерстеллар» (2014), для якого Кіп Торн і команда візуальних ефектів Double Negative розробили фізично точний GLSL-рендерер червоточини.

V. Суперформа — суперформула Гілейса на GPU

Йохан Гілейс представив суперформулу 2003 року як узагальнення кола, здатне описати надзвичайно широкий діапазон природних і математичних форм лише шістьма параметрами (m, n1, n2, n3, a, b):

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

У полярних координатах r(θ) дає радіус під кутом θ. Для 3D-суперформи два незалежні обчислення суперформули застосовуються як радіуси широти й довготи, породжуючи сфероподібну поверхню, деформовану добутком двох радіальних функцій. GPU відображає кожен піксель безпосередньо в координату (θ, φ), обчислює обидві формули, обчислює 3D-позицію на поверхні й рендерить через стандартну модель освітлення — жодної сітки, жодного геометричного шейдера, лише чиста параметрична поверхня, обчислена на піксель.

VI. Частинки GPGPU — 250 000 частинок через render-to-texture

Центральний виклик GPU-симуляції частинок — обмеження WebGL: шейдер не може одночасно читати з тієї самої текстури й писати в неї. Техніка ping-pong розв'язує це, підтримуючи дві пари цілей рендерингу (A і B). У непарних кадрах обчислювальний прохід читає з A і пише в B; у парних — читає з B і пише в A. ЦП більше не торкається даних частинок після ініціалізації.

// Псевдокод для одного кадру:
gl.bindFramebuffer(writeFBO);     // прив'язуємо ціль B
gl.useProgram(computeShader);
gl.bindTexture(readTexture);      // прив'язуємо ціль A
gl.drawArrays(TRIANGLES, 0, 6);  // повноекранний квад: запускає фізику на GPU
swap(A, B);                       // наступний кадр читає з B

// Прохід рендерингу:
gl.bindFramebuffer(null);         // задній буфер
gl.useProgram(renderShader);
gl.bindTexture(positionTexture);  // поточні позиції (щойно записані)
gl.drawArrays(POINTS, 0, N);     // кожна вершина точки вибирає свій власний тексель

Позиція кожної частинки зберігається як тексель RGBA32F (чотири 32-бітні числа з рухомою комою: позиція xyz, w для віку чи маси). Текстура 512×512 зберігає 262 144 частинки. Обчислювальний шейдер — це повноекранний квад, чий фрагментний шейдер шукає позицію частинки в gl_FragCoord, застосовує фізичне оновлення (гравітацію, шум завихрення чи атрактор) і записує нову позицію як вихідний колір. Шейдер рендерингу — це виклик малювання точок: кожен індекс вершини точки відображається на координату текселя, зчитує позицію через вибірку текстури у вершинному шейдері та виводить точковий спрайт.

Шум завихрення та бездивергентний потік

Найвізуально переконливіша поведінка частинок використовує шум завихрення (curl noise) для поля швидкостей. Для будь-якого 3D-векторного потенціалу A(p) ротор F = ∇ × A математично гарантовано є бездивергентним:

∇ · (∇ × A)  =  0   (тотожність векторного числення)

Бездивергентне поле не має джерел чи стоків — частинки ні не збиваються в щільні кульки, ні не розлітаються врізнобіч. Натомість вони формують нестисливі вихрові потоки, які разюче нагадують дим, хмари й океанські течії. Векторний потенціал A береться як 3D симплекс-шум, а ротор обчислюється аналітично (градієнт шуму диференційовний). Адвекція 250 000 частинок через шум завихрення на 60 кадрах/с не вимагає жодної участі ЦП після компіляції шейдера — GPU повністю обробляє весь конвеєр симуляції й рендерингу.