Реакція-дифузія: візерунки Тюрінга на GPU

Рівняння Грея-Скотта

Модель Грея-Скотта відстежує два хімічні види, U та V, на 2D-сітці. U постійно подається в систему; V споживається. Їхня взаємодія — реакція U + 2V → 3V:

Два параметри, f і k, контролюють усе. Крихітні зміни їхніх значень породжують радикально різні типи візерунків: плями, смуги, лабіринти, отвори, корали, відбитки пальців.

Чому CPU — неправильне місце для цього

Для сітки 512×512 це 262 144 клітинки. Кожній клітинці потрібні значення її чотирьох сусідів (лапласіан), потім два множення з плаваючою комою та два додавання. Виконання цього на JavaScript на 60 FPS вимагає ~31 мільйона операцій за секунду — досяжно, але лише ледь-ледь, і без запасу для більшої сітки чи будь-якої постобробки.

Рішення з GPU очевидне: оновлення Грея-Скотта — класична вкрай паралельна операція — кожна клітинка незалежна. Фрагментний шейдер, що читає з однієї текстури й записує в іншу, виконує те саме оновлення для всіх 262 144 клітинок одночасно.

Пінг-понг фреймбуфери

Не можна читати й писати в ту саму текстуру в одному проході рендерингу (невизначена поведінка у WebGL). Рішення — пінг-понг: підтримувати дві цілі рендерингу і чергувати, з якої читати й у яку писати кожен кадр.

// Налаштування: дві однакові цілі рендерингу
const size = 512;
const opts = {
  type: THREE.FloatType,        // float32 для точності
  format: THREE.RGBAFormat,     // R = концентрація U, G = концентрація V
  minFilter: THREE.LinearFilter,
  magFilter: THREE.LinearFilter
};
let rtA = new THREE.WebGLRenderTarget(size, size, opts);
let rtB = new THREE.WebGLRenderTarget(size, size, opts);

// Кожен кадр: обмін
function tick() {
  // Крок 1: рендеринг оновлення Грея-Скотта (читати A, писати B)
  uniforms.tState.value = rtA.texture;
  renderer.setRenderTarget(rtB);
  renderer.render(quadScene, orthoCamera);

  // Крок 2: показати результат на екрані
  displayUniforms.tState.value = rtB.texture;
  renderer.setRenderTarget(null);
  renderer.render(displayScene, orthoCamera);

  // Обмін
  [rtA, rtB] = [rtB, rtA];
  requestAnimationFrame(tick);
}

Фрагментний шейдер оновлення

GLSL-шейдер оновлення зчитує поточну текстуру стану і застосовує рівняння Грея-Скотта. Дискретний лапласіан використовує 5-точковий шаблон — класичне скінченно-різницеве наближення.

uniform sampler2D tState;   // поточна RG-текстура (R=U, G=V)
uniform vec2 uTexelSize;    // 1/ширина, 1/висота
uniform float uF;           // швидкість подачі
uniform float uK;           // швидкість знищення
uniform float uDu;          // дифузія U
uniform float uDv;          // дифузія V
uniform float uDt;          // крок часу за кадр

void main() {
  vec2 uv = vUv;
  vec2 ts = uTexelSize;

  // Вибірка поточної клітинки та 4 сусідів
  vec2 curr  = texture2D(tState, uv).rg;
  vec2 top   = texture2D(tState, uv + vec2( 0, ts.y)).rg;
  vec2 bot   = texture2D(tState, uv + vec2( 0,-ts.y)).rg;
  vec2 left  = texture2D(tState, uv + vec2(-ts.x, 0)).rg;
  vec2 right = texture2D(tState, uv + vec2( ts.x, 0)).rg;

  // 5-точковий дискретний лапласіан
  vec2 lap = top + bot + left + right - 4.0 * curr;

  float u = curr.r;
  float v = curr.g;
  float uvv = u * v * v;

  // Рівняння реакції-дифузії Грея-Скотта
  float du = uDu * lap.r - uvv + uF * (1.0 - u);
  float dv = uDv * lap.g + uvv - (uF + uK) * v;

  gl_FragColor = vec4(
    clamp(u + du * uDt, 0.0, 1.0),
    clamp(v + dv * uDt, 0.0, 1.0),
    0.0, 1.0
  );
}

Зоопарк параметрів

Найцікавіший аспект Грея-Скотта — фазова діаграма. Фіксуючи D_u = 0.21 та D_v = 0.105, ось візерунки, що виникають лише через зміну F та k:

Засівання: звідки беруться візерунки

Початковий стан має надзвичайне значення. Старт з однорідного стану не дає жодного візерунка — це стабільна рівновага. Щоб підштовхнути систему в режим формування візерунків, я засіваю невеликі «краплі» високої концентрації V у полі високого U:

// Засів початкової текстури (на CPU, завантажується один раз)
const data = new Float32Array(size * size * 4);
for (let i = 0; i < size * size; i++) {
  data[i * 4 + 0] = 1.0; // U = 1 всюди
  data[i * 4 + 1] = 0.0; // V = 0 всюди
}
// Кинути 20 випадкових насінин V
for (let s = 0; s < 20; s++) {
  const cx = Math.floor(Math.random() * size);
  const cy = Math.floor(Math.random() * size);
  for (let dx = -4; dx <= 4; dx++) {
    for (let dy = -4; dy <= 4; dy++) {
      const idx = ((cy + dy) * size + (cx + dx)) * 4;
      data[idx + 0] = 0.5; // зменшити U
      data[idx + 1] = 0.25; // додати V
    }
  }
}
const initTexture = new THREE.DataTexture(data, size, size, THREE.RGBAFormat, THREE.FloatType);

Інтерактивне малювання

Найзадовільніша функція: клацання й перетягування для введення V безпосередньо в текстуру під час виконання. Обробник події миші перетворює екранні координати в UV-простір, потім рендерить невеликий «пензлик» високого V у поточний фреймбуфер. Візерунки, які ви малюєте, органічно зростають і розвиваються, зливаючись з наявною структурою.

Це вимагало уваги до перетворень координат — квад малюється в нормалізованих координатах пристрою, але позиція миші — у CSS-пікселях. Один виклик getBoundingClientRect() та перетворення viewport-в-UV вирішили це чисто.

Що здивувало мене найбільше

Дві речі. По-перше, наскільки повільні ці візерунки в біологічних масштабах часу — симуляція виконує тисячі кроків Ейлера за секунду, тоді як реальне формування візерунків у організмі відбувається протягом годин чи днів ембріонального розвитку. Приголомшливо складна структура, яка кристалізується за секунди в симуляції, потребує повного циклу розвитку організму.

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