前言

如果你曾經仔細觀察過珊瑚的表面、河豚的花紋、或者貝殼上的條紋，你可能會好奇：這些圖案是怎麼來的？為什麼不同物種的斑紋看起來風格各異，卻又似乎遵循著某種共通的規律？

1952 年，計算機科學之父 Alan Turing 在他的論文 “The Chemical Basis of Morphogenesis” 中提出了一個驚人的假說：生物的花紋可能是由兩種化學物質的反應與擴散所形成的。這個模型後來被稱為 Reaction-Diffusion（反應擴散）系統。

更驚人的是——這個模型用程式跑出來的結果，真的和自然界的斑紋幾乎一模一樣。

這篇文章會帶你理解 Gray-Scott 模型——最流行的反應擴散模型之一，然後用 p5.js 的像素操作從零開始實作。

反應擴散的核心概念

想像一個二維平面上有兩種化學物質，我們叫它們 A 和 B。

  • A 是「食物」，它會自然地被補充（feed）進來
  • B 是「催化劑」，它會消耗 A 並繁殖自己
  • 兩者都會擴散到鄰近區域，但擴散速度不同
  • B 也會自然地衰減（kill）

用數學寫出來就是 Gray-Scott 模型的兩個偏微分方程：

∂A/∂t = Dₐ∇²A - AB² + f(1-A)
∂B/∂t = D_b∇²B + AB² - (k+f)B

其中：

  • Dₐ, D_b：A 和 B 的擴散速率（通常 Dₐ > D_b，A 擴散得比較快）
  • ∇²：拉普拉斯算子（Laplacian），代表「鄰居的平均值 – 自己的值」
  • AB²：反應項——A 和兩個 B 碰面時，A 被消耗，B 增加
  • f：feed rate，A 的補充速率
  • k：kill rate，B 的衰減速率

直覺理解

想像這樣一個生態：

  1. 食物 A 均勻地從外部灌入
  2. 催化劑 B 需要 A 才能繁殖（而且需要兩個 B 合作才能吃掉一個 A）
  3. A 擴散很快（像水一樣到處流），B 擴散很慢（像黏稠的液體）
  4. B 會自然死去

在這種情況下，B 無法均勻擴散——它在某些地方聚集、繁殖，在其他地方消亡。這種「局部的正回饋 + 全局的負回饋」就產生了各種圖案。

離散化與拉普拉斯算子

在程式中，我們用一個二維網格來近似連續的平面。每個格子存放 A 和 B 的濃度值。

拉普拉斯算子 ∇² 在離散網格上可以用卷積核來近似。最常見的是這個 3×3 的核：

| 0.05  0.2   0.05 |
| 0.2  -1.0   0.2  |
| 0.05  0.2   0.05 |

中心的權重是 -1，周圍鄰居的權重加起來也是 1，所以總和為 0。這代表「鄰居的加權平均 – 自己的值」——如果周圍濃度比自己高，拉普拉斯值為正（物質會流入），反之為負（物質會流出）。

p5.js 完整實作

let grid, next;
let dA = 1.0;
let dB = 0.5;
let feed = 0.055;
let kill = 0.062;
let w, h;


function setup() {
  pixelDensity(1);
  w = 200;
  h = 200;
  createCanvas(w, h);



// 初始化網格
  grid = [];
  next = [];
  for (let x = 0; x < w; x++) {
    grid[x] = [];
    next[x] = [];
    for (let y = 0; y < h; y++) {
      grid[x][y] = { a: 1, b: 0 };
      next[x][y] = { a: 1, b: 0 };
    }
  }



// 在中心放一些 B 作為種子
  let cx = w / 2, cy = h / 2;
  for (let x = cx - 10; x < cx + 10; x++) {
    for (let y = cy - 10; y < cy + 10; y++) {
      if (x >= 0 && x < w && y >= 0 && y < h) {
        grid[x][y].b = 1;
      }
    }
  }
}



function draw() {
  // 每幀多跑幾步以加速
  for (let step = 0; step < 10; step++) {
    simulate();
  }



// 繪製
  loadPixels();
  for (let x = 0; x < w; x++) {
    for (let y = 0; y < h; y++) {
      let idx = (x + y  w)  4;
      let a = grid[x][y].a;
      let b = grid[x][y].b;



// 著色方式：A-B 的差值
      let c = constrain(floor((a - b) * 255), 0, 255);
      pixels[idx] = c;
      pixels[idx + 1] = c;
      pixels[idx + 2] = c;
      pixels[idx + 3] = 255;
    }
  }
  updatePixels();
}



function simulate() {
  for (let x = 1; x < w - 1; x++) {
    for (let y = 1; y < h - 1; y++) {
      let a = grid[x][y].a;
      let b = grid[x][y].b;



// 計算拉普拉斯
      let laplaceA = laplacian(x, y, 'a');
      let laplaceB = laplacian(x, y, 'b');



// 反應項
      let reaction = a  b  b;



// 更新
      next[x][y].a = a + (dA  laplaceA - reaction + feed  (1 - a));
      next[x][y].b = b + (dB  laplaceB + reaction - (kill + feed)  b);



// 限制範圍
      next[x][y].a = constrain(next[x][y].a, 0, 1);
      next[x][y].b = constrain(next[x][y].b, 0, 1);
    }
  }



// 交換 grid 和 next
  let temp = grid;
  grid = next;
  next = temp;
}



function laplacian(x, y, chemical) {
  let sum = 0;
  let v = grid[x][y][chemical];



// 3x3 卷積核
  sum += grid[x-1][y-1][chemical] * 0.05;
  sum += grid[x  ][y-1][chemical] * 0.2;
  sum += grid[x+1][y-1][chemical] * 0.05;
  sum += grid[x-1][y  ][chemical] * 0.2;
  sum += v * -1.0;
  sum += grid[x+1][y  ][chemical] * 0.2;
  sum += grid[x-1][y+1][chemical] * 0.05;
  sum += grid[x  ][y+1][chemical] * 0.2;
  sum += grid[x+1][y+1][chemical] * 0.05;



return sum;
}



// 滑鼠點擊添加更多 B 種子
function mouseDragged() {
  let mx = floor(mouseX);
  let my = floor(mouseY);
  let r = 5;
  for (let dx = -r; dx <= r; dx++) {
    for (let dy = -r; dy <= r; dy++) {
      let x = mx + dx;
      let y = my + dy;
      if (x >= 0 && x < w && y >= 0 && y < h) {
        grid[x][y].b = 1;
      }
    }
  }
}

跑起來之後

一開始你只會看到中心有一小塊深色區域（B 的種子），但隨著時間推移，B 開始擴散、反應，逐漸形成各種美麗的圖案。耐心等待！通常需要幾千步迭代才能看到穩定的圖案。

feed 和 kill：圖案的密碼

Gray-Scott 模型最神奇的地方在於：只要微調 feed（f）和 kill（k）這兩個參數，就能產生截然不同的圖案。以下是一些經典的參數組合：

// 斑點（Spots）— 像河豚的斑點
let f = 0.035, k = 0.065;


// 條紋（Stripes）— 像斑馬的條紋
let f = 0.025, k = 0.056;



// 迷宮（Maze/Labyrinth）— 像珊瑚的紋路
let f = 0.029, k = 0.057;



// 脈動（Pulses）— 不斷分裂的圓點
let f = 0.014, k = 0.045;



// 蟲子（Worms）— 像蠕動的蟲
let f = 0.038, k = 0.061;



// 泡沫（Bubbles）— 像肥皂泡沫
let f = 0.012, k = 0.050;



// 漩渦（Spirals）— 旋轉的螺旋
let f = 0.018, k = 0.051;



// 混沌波動
let f = 0.026, k = 0.051;

參數探索器

加入滑桿讓你即時調整 f 和 k：

let sliderF, sliderK;


function setup() {
  // ... 之前的 setup 內容 ...



sliderF = createSlider(0, 80, 55); // f * 1000
  sliderK = createSlider(0, 80, 62); // k * 1000
  sliderF.position(220, 10);
  sliderK.position(220, 40);
}



function draw() {
  feed = sliderF.value() / 1000;
  kill = sliderK.value() / 1000;



// ... 之前的 draw 內容 ...



// 顯示參數
  fill(255);
  noStroke();
  textSize(12);
  text("f = " + nf(feed, 1, 4), 220, 75);
  text("k = " + nf(kill, 1, 4), 220, 90);
}

我強烈建議你花時間在 f-k 參數空間中探索。Robert Munafo 的 Gray-Scott 參數地圖是一個很好的參考——它把不同 f-k 值對應的圖案類型標示了出來。

效能優化

200×200 的解析度雖然能看到圖案，但有點粗糙。要提高到更大的尺寸，有幾個優化方向：

使用 Float32Array

// 用 TypedArray 取代物件陣列，效能大幅提升
let gridA, gridB, nextA, nextB;


function setup() {
  w = 400;
  h = 400;
  createCanvas(w, h);
  pixelDensity(1);



gridA = new Float32Array(w * h).fill(1);
  gridB = new Float32Array(w * h).fill(0);
  nextA = new Float32Array(w * h);
  nextB = new Float32Array(w * h);



// 初始種子
  let cx = w / 2, cy = h / 2;
  for (let dx = -10; dx < 10; dx++) {
    for (let dy = -10; dy < 10; dy++) {
      let idx = (cx + dx) + (cy + dy) * w;
      gridB[idx] = 1;
    }
  }
}



function simulate() {
  for (let y = 1; y < h - 1; y++) {
    for (let x = 1; x < w - 1; x++) {
      let idx = x + y * w;
      let a = gridA[idx];
      let b = gridB[idx];



// 拉普拉斯（內聯以減少函數呼叫開銷）
      let lapA =
        gridA[idx - w - 1] * 0.05 +
        gridA[idx - w]     * 0.2  +
        gridA[idx - w + 1] * 0.05 +
        gridA[idx - 1]     * 0.2  +
        a                  * -1.0 +
        gridA[idx + 1]     * 0.2  +
        gridA[idx + w - 1] * 0.05 +
        gridA[idx + w]     * 0.2  +
        gridA[idx + w + 1] * 0.05;



let lapB =
        gridB[idx - w - 1] * 0.05 +
        gridB[idx - w]     * 0.2  +
        gridB[idx - w + 1] * 0.05 +
        gridB[idx - 1]     * 0.2  +
        b                  * -1.0 +
        gridB[idx + 1]     * 0.2  +
        gridB[idx + w - 1] * 0.05 +
        gridB[idx + w]     * 0.2  +
        gridB[idx + w + 1] * 0.05;



let reaction = a  b  b;



nextA[idx] = a + dA  lapA - reaction + feed  (1 - a);
      nextB[idx] = b + dB  lapB + reaction - (kill + feed)  b;
    }
  }



// 交換
  let tempA = gridA; gridA = nextA; nextA = tempA;
  let tempB = gridB; gridB = nextB; nextB = tempB;
}

使用 WebGL (GLSL) 加速

最終極的方案是把整個模擬搬到 GPU 上。用兩張浮點紋理做 ping-pong buffer，在 fragment shader 中進行迭代。這樣即使 1024×1024 的解析度也能跑到 60fps：

precision highp float;


uniform sampler2D u_texture; // 上一幀的狀態
uniform vec2 u_resolution;
uniform float u_feed;
uniform float u_kill;



void main() {
    vec2 uv = gl_FragCoord.xy / u_resolution;
    vec2 texel = 1.0 / u_resolution;



// 讀取當前值
    vec4 val = texture2D(u_texture, uv);
    float a = val.r;
    float b = val.g;



// 拉普拉斯（讀取鄰居像素）
    float lapA = 0.0, lapB = 0.0;



for (int dy = -1; dy <= 1; dy++) {
        for (int dx = -1; dx <= 1; dx++) {
            vec2 offset = vec2(float(dx), float(dy)) * texel;
            vec4 neighbor = texture2D(u_texture, uv + offset);



float weight;
            if (dx == 0 && dy == 0) weight = -1.0;
            else if (dx == 0 || dy == 0) weight = 0.2;
            else weight = 0.05;



lapA += neighbor.r * weight;
            lapB += neighbor.g * weight;
        }
    }



float reaction = a  b  b;
    float dA = 1.0, dB = 0.5;



float newA = a + dA  lapA - reaction + u_feed  (1.0 - a);
    float newB = b + dB  lapB + reaction - (u_kill + u_feed)  b;



gl_FragColor = vec4(clamp(newA, 0.0, 1.0), clamp(newB, 0.0, 1.0), 0.0, 1.0);
}

著色：讓圖案更好看

黑白的反應擴散已經很美了，但加上色彩能讓它更驚艷：

// 在 draw() 的著色部分
let a = gridA[idx];
let b = gridB[idx];


// 方式一：冷暖色調
let r = constrain(floor((1 - b)  200 + b  80), 0, 255);
let g = constrain(floor((1 - b)  200 + b  40), 0, 255);
let blue = constrain(floor(b * 255), 0, 255);



// 方式二：火焰色
let t = b;
let r2 = constrain(floor(t  3  255), 0, 255);
let g2 = constrain(floor(t  t  255), 0, 255);
let b2 = constrain(floor(pow(t, 3) * 255), 0, 255);



// 方式三：生物感
let base = (a - b);
let r3 = constrain(floor(base * 180 + 50), 0, 255);
let g3 = constrain(floor(base * 200 + 40), 0, 255);
let b3 = constrain(floor(base * 100 + 70), 0, 255);

與自然的連結

Alan Turing 在 1952 年提出反應擴散假說時，科學界半信半疑。但近年來的研究越來越支持他的理論：

  • 斑馬魚的條紋確實由兩種色素細胞的擴散和相互作用形成
  • 貝殼上的花紋可以用一維反應擴散模型精確模擬（因為貝殼是從邊緣一行一行地長出來的）
  • 珊瑚的表面紋理與 Gray-Scott 模型的迷宮模式驚人地相似
  • 手指指紋的形成也被認為涉及類似的機制

這是一個美麗的案例：數學模型不只是「長得像」自然，而是真的揭示了自然的運作機制。

進階玩法

非均勻參數

讓 feed 和 kill 在空間上變化，可以在同一張圖中看到不同類型的圖案：

// 在 simulate() 中，根據位置調整 f 和 k
let localFeed = feed + (x / w - 0.5) * 0.02;
let localKill = kill + (y / h - 0.5) * 0.01;

多種子隨機初始化

// 隨機散佈多個 B 種子
for (let i = 0; i < 20; i++) {
  let sx = floor(random(20, w - 20));
  let sy = floor(random(20, h - 20));
  for (let dx = -5; dx <= 5; dx++) {
    for (let dy = -5; dy <= 5; dy++) {
      let idx = (sx + dx) + (sy + dy) * w;
      gridB[idx] = 1;
    }
  }
}

小結

反應擴散系統是一個絕佳的例子，展示了簡單的局部規則如何產生全局的複雜圖案。兩種化學物質、兩個擴散速率、兩個控制參數——就這麼少的要素，卻能模擬出自然界中千變萬化的斑紋和圖案。

Turing 在 1952 年的這篇論文，不僅開創了數學生物學的新領域，也給了我們程式藝術家一個強大的工具。每次看到反應擴散的模擬結果，我都會感歎：自然界最深層的美學，竟然可以用幾行方程式來捕捉。

延伸閱讀

  • Karl Sims 的 “Reaction-Diffusion Tutorial”——最經典的入門教學
  • Robert Munafo 的 Gray-Scott 參數地圖 (mrob.com/pub/comp/xmorphia)
  • Turing 的原始論文 “The Chemical Basis of Morphogenesis” (1952)
  • Daniel Shiffman Coding Challenge #13: Reaction Diffusion Algorithm
  • NVIDIA 的 CUDA 實作範例：用 GPU 跑超高解析度的反應擴散