Multi Scale Turing Pattern

目次

1. はじめに

2. Multi Scale Turing Patternとは

3. アルゴリズムの概要

4. 実装

5. まとめ

6. 参考

はじめに

こんにちは、NUMです。 今回はMulti Scale Turing Patternという少し変わったTuring Patternについて紹介します。

Multi Scale Turing Patternとは

Turing Patternは反応拡散方程式とも呼ばれ、活性因子と抑制因子の相互作用によって複雑な空間パターンを生み出します。Turing Patternについては以前記事を書いているのでぜひ読んでみて下さい。

Turing Pattern

元々のTuring Patternはこのように模様が1層しかないですが、Multi Scale Turing Patternは層が複数存在するように見え、複雑かつ生物らしいTuring Patternです。生物というのは細胞が多層構造で構成されているので、それらしく見えるのは不思議なことではないと思います。

Multi Scale Turing Pattern

アルゴリズムの概要

今回の反応拡散方程式を数値的に解くのではなく、Box Blarという拡散を近似するアルゴリズムを用いて、各スケールごとの活性因子と抑制因子の拡散を計算し、局所で「どのスケールが最も安定か」を判定してグリッドを更新する方式をとっています。

処理の流れ

1. 3次元配列のactivator, inhibitor, variationを用意 activator：活性因子の場 inhibitor：抑制因子の場 variation：活性因子と抑制因子の差分値を格納する場
   

   

   
   

2. 各スケールで活性因子と抑制因子の拡散をBox Blurで拡散を計算 Turing Patternの拡散項に相当
   

   

   
   

3. 各スケールで全要素の活性因子と抑制因子の差を計算
   

   

4. 各ピクセルで variation が最小のスケール bestScale を選ぶ variationが最小ということは安定した模様になっている
   

   

5. bestScale の活性因子、抑制因子の大小に応じてgrid[x][y]の値を増加、減少させる Turing Patternの反応項に相当 grid は2次元配列なので各スケールの相互作用の結果が1面に反映されているイメージ
   

   

   
   

6. 最後に grid 全体を-1~1に正規化
   

実装

int SCALES = 4; // スケール数

float[][]   grid; // メイングリッド(時刻tにおける物質の量) [x][y]
float[][][] activator, inhibitor, variation; // 各スケール用の活性因子/抑制因子 [scale][x][y]
int[] actRadius, inhRadius; // 活性因子/抑制因子の拡散半径
float[] smallAmount; // 物質の更新量

color[] scaleColor; // カラー出力時のスケール色
color[][] colGrid;  // 選択したスケール色を保持

void setup() {
  size(800, 800);
  pixelDensity(displayDensity());
  background(0);
  
  initialize();
  for (int i = 0; i < 200; i++) {stepSimulation();}
  renderToScreen();
}

/**
 * 配列の初期化、メイングリッドの初期値代入
 */
void initialize() {
  grid = new float[width][height];

  activator = new float[SCALES][width][height];
  inhibitor = new float[SCALES][width][height];
  variation = new float[SCALES][width][height];

  actRadius = new int[] {4, 20, 100, 500};
  inhRadius = new int[] {8, 40, 200, 1000};
  smallAmount = new float[] {0.05, 0.04, 0.03, 0.02};

  scaleColor = new color[] {color(0, 30, 0), color(20, 80, 20), color(100, 180, 100), color(210, 250, 210)};
  colGrid = new color[width][height];

  // メイングリッドを-1~1の範囲で初期化
  for (int x = 0; x < width; x++) {
    for (int y = 0; y < height; y++) {
      grid[x][y] = random(-1, 1);
    }
  }
}

/**
 * ①拡散：活性因子と抑制因子のボックスブラーで計算
 * ②活性因子から抑制因子の差分を算出
 * ③差分から一番安定したスケールを選定
 * ④反応：活性因子と抑制因子を比較して物質を増減
 * ⑤メイングリッドを-1~1に正規化
 */
void stepSimulation() {
  // 各スケールで活性因子/抑制因子を計算
  for (int s = 0; s < SCALES; s++) {
    boxBlur(grid, activator[s], actRadius[s]);
    boxBlur(grid, inhibitor[s], inhRadius[s]);
  }

  // 後続の処理でvariationが最小のスケールを取得するための事前準備
  for (int s = 0; s < SCALES; s++) {
    for (int x = 0; x < width; x++) {
      for (int y = 0; y < height; y++) {
        variation[s][x][y] = abs(activator[s][x][y] - inhibitor[s][x][y]);
      }
    }
  }

  // 各ピクセルを更新：variationが最小のスケールを選び、そのスケールの活性因子/抑制因子によって物質を増減
  for (int x = 0; x < width; x++) {
    for (int y = 0; y < height; y++) {
      int bestScale = 0;
      float bestVar = variation[0][x][y];

      // 安定したパターンの選定
      // 各スケールの同じ位置にあるvariationをチェックし、最小のvariationを持つスケールを取得
      for (int s = 1; s < SCALES; s++) {
        if (variation[s][x][y] < bestVar) {
          bestVar = variation[s][x][y];
          bestScale = s;
        }
      }

      // 反応項：活性因子が抑制因子より多いと物質を増加させる
      if (activator[bestScale][x][y] > inhibitor[bestScale][x][y]) {
        grid[x][y] += smallAmount[bestScale];
      } else {
        grid[x][y] -= smallAmount[bestScale];
      }
      
      // 安定しているスケールの可視化をするため採用する色をスケール毎に変える
      colGrid[x][y] = scaleColor[bestScale];
    }
  }

  //メイングリッドを-1~1の範囲に正規化
  normalizeGrid();
}

/**
 * 画像の二次元配列に対してボックスブラー（移動平均フィルタ）を適用する
 *
 * ブラー処理は2段階で行い、
 * 水平方向の移動平均を計算し、一時配列tmpに保存
 * 垂直方向の移動平均を計算し、結果をdstに保存
 *
 * @param src    入力配列(メイングリッド)
 * @param dst    出力配列(活性因子/抑制因子)
 * @param radius ブラーの半径
 */
void boxBlur(float[][] src, float[][] dst, int radius) {
  float[][] tmp = new float[width][height]; // 平均値を格納
  int window = radius * 2 + 1;

  // 水平セルの処理
  for (int y = 0; y < height; y++) {
    float sum = 0;

    // 水平セルの座標(0,y)の平均
    for (int i = -radius; i <= radius; i++) {
      int xi = constrain(i, 0, width - 1);
      sum += src[xi][y];
    }
    tmp[0][y] = sum / window;

    // 合計値は差分のみ計算して取得し、平均をtmpに格納
    // x=0: [-2, -1, 0, 1, 2]
    // x=1: [-1,  0, 1, 2, 3]
    // x=0: [-2, -1, 0, 1, 2]からaddIdx = 3を追加、subIdx = -2を引くことで差分更新
    for (int x = 1; x < width; x++) {
      int addIdx = x + radius;
      int subIdx = x - radius - 1;
      addIdx = constrain(addIdx, 0, width - 1);
      subIdx = constrain(subIdx, 0, width - 1);
      sum += src[addIdx][y] - src[subIdx][y];
      tmp[x][y] = sum / window;
    }
  }

  // 垂直セルの処理
  for (int x = 0; x < width; x++) {
    float sum = 0;

    // 垂直セルの座標(0,y)の平均
    for (int j =- radius; j <= radius; j++) {
      int yj = constrain(j, 0, height - 1);
      sum += tmp[x][yj];
    }
    dst[x][0] = sum / window;

    // 垂直セルの平均値をdstに格納(合計値の計算方法は水平と同じ)
    for (int y = 1; y < height; y++) {
      int addIdx = y + radius;
      int subIdx = y - radius - 1;
      addIdx = constrain(addIdx, 0, height - 1);
      subIdx = constrain(subIdx, 0, height - 1);
      sum += tmp[x][addIdx] - tmp[x][subIdx];
      dst[x][y] = sum / window;
    }
  }
}

/**
 * メイングリッドを-1~1の範囲に正規化
 */
void normalizeGrid() {
  float minv = Float.MAX_VALUE;
  float maxv = -Float.MAX_VALUE;

  // メイングリッドの最小値、最大値を抽出
  for (int x = 0; x < width; x++) {
    for (int y = 0; y < height; y++) {
      float v = grid[x][y];
      if (v < minv) minv = v;
      if (v > maxv) maxv = v;
    }
  }

  // 範囲の取得
  float range = maxv - minv;
  if (range == 0) return;

  // ① 各要素に最小値を引いて0~rangeにシフト
  // ② 範囲で割って0~1の範囲にスケール
  // ③ 2を掛けて0~2に拡げる
  // ④ -1を足して-1~1にシフト
  for (int x = 0; x < width; x++) {
    for (int y = 0; y < height; y++) {
      grid[x][y] = ((grid[x][y] - minv) / range) * 2.0 - 1.0;
    }
  }
}

/**
 * メイングリッドの可視化
 */
void renderToScreen() {
  noStroke();
  for (int x = 0; x < width; x++) {
    for (int y = 0; y < height; y++) {
      float v = (grid[x][y] + 1.0) * 0.5; // -1~1を0~1に変換
      float r = red(colGrid[x][y]) * v;
      float g = green(colGrid[x][y]) * v;
      float b = blue(colGrid[x][y]) * v;
      color c = color(constrain(r, 0, 255), constrain(g, 0, 255), constrain(b, 0, 255));
      fill(c);
      rect(x, y, 1, 1);
    }
  }
}

初期化

/**
 * 配列の初期化、メイングリッドの初期値代入
 */
void initialize() {
  grid = new float[width][height];

  activator = new float[SCALES][width][height];
  inhibitor = new float[SCALES][width][height];
  variation = new float[SCALES][width][height];

  actRadius = new int[] {4, 20, 100, 500};
  inhRadius = new int[] {8, 40, 200, 1000};
  smallAmount = new float[] {0.05, 0.04, 0.03, 0.02};

  scaleColor = new color[] {color(0, 30, 0), color(20, 80, 20), color(100, 180, 100), color(210, 250, 210)};
  colGrid = new color[width][height];

  // メイングリッドを-1~1の範囲で初期化
  for (int x = 0; x < width; x++) {
    for (int y = 0; y < height; y++) {
      grid[x][y] = random(-1, 1);
    }
  }
}

initialize関数では各配列のの初期化、値の代入を行います。 grid配列には-1~1の範囲でランダムな値が代入され、後続のstepSimulation関数で値が増減されていきます。

Box Blur

/**
 * 画像の二次元配列に対してボックスブラー（移動平均フィルタ）を適用する
 *
 * ブラー処理は2段階で行い、
 * 水平方向の移動平均を計算し、一時配列tmpに保存
 * 垂直方向の移動平均を計算し、結果をdstに保存
 *
 * @param src    入力配列(メイングリッド)
 * @param dst    出力配列(活性因子/抑制因子)
 * @param radius ブラーの半径
 */
void boxBlur(float[][] src, float[][] dst, int radius) {
  float[][] tmp = new float[width][height]; // 平均値を格納
  int window = radius * 2 + 1;

  // 水平セルの処理
  for (int y = 0; y < height; y++) {
    float sum = 0;

    // 水平セルの座標(0,y)の平均
    for (int i = -radius; i <= radius; i++) {
      int xi = constrain(i, 0, width - 1);
      sum += src[xi][y];
    }
    tmp[0][y] = sum / window;

    // 合計値は差分のみ計算して取得し、平均をtmpに格納
    // x=0: [-2, -1, 0, 1, 2]
    // x=1: [-1,  0, 1, 2, 3]
    // x=0: [-2, -1, 0, 1, 2]からaddIdx = 3を追加、subIdx = -2を引くことで差分更新
    for (int x = 1; x < width; x++) {
      int addIdx = x + radius;
      int subIdx = x - radius - 1;
      addIdx = constrain(addIdx, 0, width - 1);
      subIdx = constrain(subIdx, 0, width - 1);
      sum += src[addIdx][y] - src[subIdx][y];
      tmp[x][y] = sum / window;
    }
  }

  // 垂直セルの処理
  for (int x = 0; x < width; x++) {
    float sum = 0;

    // 垂直セルの座標(0,y)の平均
    for (int j =- radius; j <= radius; j++) {
      int yj = constrain(j, 0, height - 1);
      sum += tmp[x][yj];
    }
    dst[x][0] = sum / window;

    // 垂直セルの平均値をdstに格納(合計値の計算方法は水平と同じ)
    for (int y = 1; y < height; y++) {
      int addIdx = y + radius;
      int subIdx = y - radius - 1;
      addIdx = constrain(addIdx, 0, height - 1);
      subIdx = constrain(subIdx, 0, height - 1);
      sum += tmp[x][addIdx] - tmp[x][subIdx];
      dst[x][y] = sum / window;
    }
  }
}

Box Blurとは「近傍ピクセルの平均を取る」アルゴリズムですが、実はこれが拡散の離散的近似になっています。この処理はstepSimulation関数で呼び出されています。 この処理は文章で説明するより図の方が理解しやすいと思います。

拡散、反応、最小スケールの選定

/**
 * ①拡散：活性因子と抑制因子のボックスブラーで計算
 * ②活性因子から抑制因子の差分を算出
 * ③差分から一番安定したスケールを選定
 * ④反応：活性因子と抑制因子を比較して物質を増減
 * ⑤メイングリッドを-1~1に正規化
 */
void stepSimulation() {
  // 各スケールで活性因子/抑制因子を計算
  for (int s = 0; s < SCALES; s++) {
    boxBlur(grid, activator[s], actRadius[s]);
    boxBlur(grid, inhibitor[s], inhRadius[s]);
  }

  // 後続の処理でvariationが最小のスケールを取得するための事前準備
  for (int s = 0; s < SCALES; s++) {
    for (int x = 0; x < width; x++) {
      for (int y = 0; y < height; y++) {
        variation[s][x][y] = abs(activator[s][x][y] - inhibitor[s][x][y]);
      }
    }
  }

  // 各ピクセルを更新：variationが最小のスケールを選び、そのスケールの活性因子/抑制因子によって物質を増減
  for (int x = 0; x < width; x++) {
    for (int y = 0; y < height; y++) {
      int bestScale = 0;
      float bestVar = variation[0][x][y];

      // 安定したパターンの選定
      // 各スケールの同じ位置にあるvariationをチェックし、最小のvariationを持つスケールを取得
      for (int s = 1; s < SCALES; s++) {
        if (variation[s][x][y] < bestVar) {
          bestVar = variation[s][x][y];
          bestScale = s;
        }
      }

      // 反応項：活性因子が抑制因子より多いと物質を増加させる
      if (activator[bestScale][x][y] > inhibitor[bestScale][x][y]) {
        grid[x][y] += smallAmount[bestScale];
      } else {
        grid[x][y] -= smallAmount[bestScale];
      }
      
      // 安定しているスケールの可視化をするため採用する色をスケール毎に変える
      colGrid[x][y] = scaleColor[bestScale];
    }
  }

  //メイングリッドを-1~1の範囲に正規化
  normalizeGrid();
}

stepSimulation関数ではTuring Patternの重要な仕組みである「拡散と反応」、そしてMulti Scale Turing Pattern独自の仕組みとして最も安定したパターンの選定が行われています。 最後にgrid配列の要素を-1~1に正規化しています。

各スケールで活性因子と抑制因子の拡散を計算し、activator、inhibitor配列に格納しています。 この拡散の計算はboxBlurという高速なアルゴリズム使用して計算しています。詳細は次で解説します。

for (int s = 0; s < SCALES; s++) {
  boxBlur(grid, activator[s], actRadius[s]);
  boxBlur(grid, inhibitor[s], inhRadius[s]);
}

活性因子と抑制因子の差分を対応するスケールの要素に格納しています。 これは「活性と抑制がどれだけ乖離しているか」を表していて、乖離しているほど綺麗な模様ができにくいです。 variation[s][x][y] = abs(activator[s][x][y] - inhibitor[s][x][y]);

variationの値は同じ要素の各スケールと比較され、最小の要素を持つスケールを選んでいます。これは活性と抑制が近い（差が小さい）スケールを「その場所の安定なスケール」と見なし、そのスケールでgrid配列を更新するためです。

for (int s = 0; s < SCALES; s++) {
  for (int x = 0; x < width; x++) {
    for (int y = 0; y < height; y++) {
      variation[s][x][y] = abs(activator[s][x][y] - inhibitor[s][x][y]);
    }
  }
}

通常のTuring Patternにおける反応項に相当する処理です。 一番安定しているスケールの活性因子が抑制因子より多ければ値を増加、少なければ値を減少させます。更新量はスケール毎に変えています。

if (activator[bestScale][x][y] > inhibitor[bestScale][x][y]) {
  grid[x][y] += smallAmount[bestScale];
} else {
  grid[x][y] -= smallAmount[bestScale];
}

最後にgrid配列の値をnormalizeGrid()で-1~1に正規化します。 これは正規化することで反応項によって値が無限に増加もしくは減少してしまうことを防ぎます。

正規化

/**
 * メイングリッドを-1~1の範囲に正規化
 */
void normalizeGrid() {
  float minv = Float.MAX_VALUE;
  float maxv = -Float.MAX_VALUE;

  // メイングリッドの最小値、最大値を抽出
  for (int x = 0; x < width; x++) {
    for (int y = 0; y < height; y++) {
      float v = grid[x][y];
      if (v < minv) minv = v;
      if (v > maxv) maxv = v;
    }
  }

  // 範囲の取得
  float range = maxv - minv;
  if (range == 0) return;

  // ① 各要素に最小値を引いて0~rangeにシフト
  // ② 範囲で割って0~1の範囲にスケール
  // ③ 2を掛けて0~2に拡げる
  // ④ -1を足して-1~1にシフト
  for (int x = 0; x < width; x++) {
    for (int y = 0; y < height; y++) {
      grid[x][y] = ((grid[x][y] - minv) / range) * 2.0 - 1.0;
    }
  }
}

normalizeGrid関数はメイングリッドの要素を-1~1の範囲に正規化します。 まず全要素から最小値と最大値を取得し、最小値と最大値の範囲を計算します。 範囲を計算したら以下手順で正規化を行なっていきます。

例：grid = {2, 4, 6} のとき

• minv=2, maxv=6, range=4

• v - minv → {0, 2, 4}（0〜range）

• /range → {0, 0.5, 1.0}（0〜1）

• ×2 -1 → {-1, 0, 1}（-1〜1）
  

まとめ

今回はMulti Scale Turing Patternについて解説しました。難解そうに見えますが、通常のTuring Patternとそれほど違いはなく、各スケールの値を比較して最も安定したスケールの更新量を使用しているだけです。でも出来上がる模様は複雑で、より生命らしさに近づいているような感じですよね。

この作品はMulti Scale Turing Patternを使用して制作しました。 作品名は「Ouroboros - “死と再生”」と言い、Ouroboros（ウロボロス」をテーマとしています。 Ouroborosは自分の尾を噛んで円環となった蛇または竜を描いた古代のシンボルのことで、「循環性」「永続性」「始原性」「無限性」「完全性」といった様々意味を持ちます。

Ouroboros - "死と再生"

Multi Scale Turing Patternの模様を見たときに真っ先にこのコンセプトが合うだろうと思いました。より生命らしさを持つ模様がコードという無限の可能性をもつツールによって生み出されるというのがOuroborosとマッチしていると思いました。 今後もOuroborosシリーズの作品を作っていこうと思うのでぜひ見て下さい！ 最後まで読んでいただきありがとうございました！

参考

https://rreusser.github.io/multiscale-turing-pattern-gallery/

https://softologyblog.wordpress.com/2011/07/05/multi-scale-turing-patterns/

https://qiita.com/UWATechnology/items/c63032f0770bdaff2e24