こんにちは、NUMです。
少しづつ暖かくなってきて過ごしやすい季節になりましたね。
今回は動物の模様を作り出すチューリングパターンについて自分なりに調べたことを書いていこうと思います。
チューリングパターンはProcessingを始めた時からずっと実装してみたいと思っていました。ついにこの時が来たって感じです！
目次
チューリングパターンとは
チューリング・パターンは、数学者のアラン・チューリングによって考案された反応拡散方程式を用いたパターン形成の現象です。
例としてよく魚の模様が挙げられます。以下のようなウネウネした模様です。この有機的な模様は、ある化学物質が相互に作用することで作られることが分かっています。
反応拡散方程式とは、以下拡散方程式と反応項を組み合わせた式となります。
拡散方程式：物質の拡散や、熱の拡散に関する式
反応項：化学反応に関する項
ちなみに、チューリングのことが知りたくて、「イミテーションゲーム」という映画を見てみましたが、数学が戦争終結につながるとは思いもしませんでした。
拡散方程式
個人的に反応拡散方程式は、拡散方程式から理解していくと頭に入りやすいです。
拡散方程式は、物質の拡散や熱の拡散に関する式です。
水に絵の具を垂らした時に色が広がっていくイメージです。
拡散方程式（一次元）は上記の数式です。
u（物質の濃度を表す）が多い場合に減少して、少ない場合に増加して均一に広がっていくという内容です。
拡散方程式は、関数uをxに関して2階微分した式にD：拡散係数を掛け合わせています。
この式をプログラム化したいですが、コンピュータは微分ができないので四則演算の式に変換する必要があります。そこで「差分法」という手法で微分に近似した値を計算します。
まず微分について簡単に説明します。
微分の目的は、関数の微小な変化を分析することです。
「微小な変化」とは「傾き」のことです。傾きは「yの変化量/xの変化量」で表現できましたね。
hがxの変化量、f(a+h)-f(a)がyの変化量です。
微分は関数の微小な変化を見るために、h(xの変化量)を限りなく0に近づけます。
limΔ→0の部分がそれにあたります。
これは接線の傾きを求めることと同義になります。
微分の詳細は予備ノリ先生が分かりやすく解説しているので参照してください笑
次はコンピュータに微分を解かせるための「差分法」について説明します。
差分法は、関数のある区間の計算結果を刻み幅(Δ)で割ることで微分に近似することができる手法です。
以下が差分法の式です。Δは微小な刻み幅を表し、xの点とその近くの点 x+Δにおける関数の値を用いて微分を近似します。
それでは拡散方程式を解いていきます。
今回は差分法の1種である中心差分法という手法で解いていきます。
中心差分法は関数値を中心の点の前後の点（ui+1,ui-1）で近似し、その差を刻み幅で割ることで微分を求めます。
1. ui+1,ui-1をテイラー展開し、2階微分の項だけを取得します。
2. ui+1,ui-1を足し合わせる。O(Δx^4)は高次の項は誤差の影響が少ないので無視するという意味です。
3. 2階微分の式を整理して中心差分法の式に変換する。
これで微分を四則演算で表現できました。それではプログラムを実装していきます。
1次元拡散方程式：
int n = 200;
int stepSize;
int denominator;
float[] u = new float[n];
float[] nextU = new float[n];
float D = 0.1;
void setup() {
size(400, 400);
background(255);
pixelDensity(displayDensity());
stepSize = width/n; // 刻み幅
denominator = stepSize*stepSize; // 分母
// uの濃度をランダムに設定
for (int i = 0; i < n; i++) {
u[i] = random(0, 1);
}
}
void draw() {
background(255);
// 中心差分法は両端の値を使用する。配列の要素外参照エラーを防ぐためiを1から開始
for (int i = 1; i < n-1; i++) {
// x座標ラプラシアン
float laplacianX = (u[i+1] - 2*u[i] + u[i-1]) / denominator;
// 拡散方程式
nextU[i] = u[i] + D * (laplacianX);
}
// 描画
for (int i = 0; i < n; i++) {
// 高さにuを掛け合わせることで拡散の様相を可視化
point(i * stepSize, width/2 + (width/2 * nextU[i]));
}
// 配列uに計算後の配列nextUを上書きすることで、配列を使い回せる
System.arraycopy(nextU, 0, u, 0, n);
}
以下プログラムが特に重要な箇所です。
// x座標ラプラシアン
float laplacianX = (u[i+1] - 2*u[i] + u[i-1]) / denominator;
// 拡散方程式
nextU[i] = u[i] + D * (laplacianX);
laplacianX は、中心点の前後の点との差を取り、分母 denominator で割っています。これによって中心点の周囲の濃度変化を表すラプラシアン(※1)が求められます。
Dは拡散係数で、物質の拡散速度を示します。この値を laplacianX に掛けることで、その位置の濃度変化がどれだけ拡散によって影響を受けるかを計算します。
最終的に、u[i] に D * laplacianX を加えることで、その位置の次の時刻の濃度 nextU[i] を求めています。
※1 ラプラシアンとは、多変数関数の偏微分の和です。
偏微分とは多変数関数において、一つの変数に着目し、残りの変数は定数とみなして微分をすることです。
1次元拡散方程式だと変数はxのみなので、以下で言うとxの式だけになります。
記号は∆,∇^2で表されたりします。
2次元拡散方程式：
int n = 200;
int stepSize;
int denominator;
float[][] u = new float[n][n];
float[][] nextU = new float[n][n];
float D = 0.1;
void setup() {
size(400, 400);
pixelDensity(displayDensity());
frameRate(60);
stepSize = width/n; // 刻み幅
denominator = stepSize*stepSize; // 分母
// uの濃度をランダムに設定
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
u[i][j] = random(0, 1);
}
}
}
void draw() {
// 数値解法
for (int i = 1; i < n-1; i++) {
for (int j = 1; j < n-1; j++) {
// x座標ラプラシアン
float laplacianX = (u[i+1][j] - 2*u[i][j] + u[i-1][j]) / denominator;
// y座標ラプラシアン
float laplacianY = (u[i][j+1] - 2*u[i][j] + u[i][j-1]) / denominator;
// 拡散方程式の数値解法
nextU[i][j] = u[i][j] + D * (laplacianX + laplacianY);
}
}
// 結果の描画
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
fill(255, nextU[i][j]*255, 0);
noStroke();
rect(i*stepSize, j*stepSize, stepSize, stepSize);
}
}
// 配列uに配列nextUを上書き
System.arraycopy(nextU, 0, u, 0, n);
}
2次元でも考え方は同じです。今回はyに関する微分も使用します。
yの場合、中心差分法は近傍の上下の値を使用して計算します。
// x座標ラプラシアン
float laplacianX = (u[i+1][j] - 2*u[i][j] + u[i-1][j]) / denominator;
// y座標ラプラシアン
float laplacianY = (u[i][j+1] - 2*u[i][j] + u[i][j-1]) / denominator;
// 拡散方程式の数値解法
nextU[i][j] = u[i][j] + D * (laplacianX + laplacianY);
反応拡散方程式
それでは本題の反応拡散方程式について説明します。
今回はGray-Scottモデルを使おうと思います。これは物質U,Vの濃度を表す変数u,vの変化をモデル化したものです。
uとvはお互いに反応し、拡散する。それが時間経過に伴って、u,vの濃度に差が出てきます。これを可視化すると様々なパターンが出来上がるというイメージです。
式は一見難解に見えますが、先ほどの拡散方程式を学習したので、一つ一つ分解すれば単純な式に見えますね。
+f(1-u)：流入項
-v(f+k)：流出項
uv^2：反応項
f、k：定数
Du、Dv：拡散係数
∆ ：拡散を表現（ラプラシアン）
式が理解できたら実装していきましょう。以下概要ですが、作りは拡散方程式とほぼ同じです。
u,v,ラプラシアン用の配列を準備
特定の区間にvを設定（色を変化させたい場所）
u,vを使用してラプラシアンの計算
反応拡散方程式の計算
描画
チューリングパターンプログラム：
float[][] u, v, nextU, nextV;
float[][] laplacianXu, laplacianYu, laplacianXv, laplacianYv;
// チューリングパターン定数
float Du = 0.2; // u拡散係数
float Dv = 0.1; // v拡散係数
float f = 0.046; // 餌
float k = 0.063; // キルレート
void setup() {
size(400, 400);
pixelDensity(displayDensity());
// 成分u,v,ラプラシアン配列
u = new float[width][height];
v = new float[width][height];
nextU = new float[width][height];
nextV = new float[width][height];
laplacianXu = new float[width][height];
laplacianYu = new float[width][height];
laplacianXv = new float[width][height];
laplacianYv = new float[width][height];
// uを全要素に設定
for (int i = 0; i < width; i++) {
for (int j = 0; j < height; j++) {
u[i][j] = 1;
nextU[i][j] = 1;
// 特定の区間にvを設定(ここの条件を変えれば任意の形を作れる)
if ((i > width/2 - 20 & i < width/2 + 20) & (j > height/2 - 20 & j < height/2 + 20)) {
v[i][j] = 1;
nextV[i][j] = 1;
} else {
v[i][j] = 0;
nextV[i][j] = 0;
}
}
}
}
void draw() {
// ラプラシアンの計算
for (int i = 1; i < width-1; i ++) {
for (int j = 1; j < height-1; j ++) {
laplacianXu[i][j] = u[i+1][j] - 2*u[i][j] + u[i-1][j];
laplacianYu[i][j] = u[i][j+1] - 2*u[i][j] + u[i][j-1];
laplacianXv[i][j] = v[i+1][j] - 2*v[i][j] + v[i-1][j];
laplacianYv[i][j] = v[i][j+1] - 2*v[i][j] + v[i][j-1];
}
}
for (int i = 1; i < width-1; i++) {
for (int j = 1; j < height-1; j++) {
float laplacianU = laplacianXu[i][j] + laplacianYu[i][j];
float laplacianV = laplacianXv[i][j] + laplacianYv[i][j];
// パフォーマンス向上のため同じ計算は先にしておく
float uvv = u[i][j]*v[i][j]*v[i][j];
// 反応拡散方程式
u[i][j] += (Du * laplacianU) - uvv + f*(1-u[i][j]);
v[i][j] += (Dv * laplacianV) + uvv - (k+f)*v[i][j];
}
}
noStroke();
for (int i = 0; i < width; i++) {
for (int j = 0; j < height; j++) {
fill(u[i][j]*255);
rect(i, j, 1, 1);
}
}
// 次の計算に計算結果を使い回す
System.arraycopy(nextU, 0, u, 0, width);
System.arraycopy(nextV, 0, v, 0, width);
}
実行すると徐々に模様が広がっていき最終的にこんな模様が出来上がります。
f,kの値を変えると模様が変わりますが、良い感じの数値にしないと消えたりするので難しいです。
まとめ
今回はチューリングパターンを拡散方程式の観点から理解し、実装していくという記事を書きました。
難しい数式が出てきたので、Processing初心者の方には少々難解な内容だと思いますが
仕組みを理解することで、作品への応用ができるようになるのでオススメです！
僕は元々生物学を学んでいましたが、数学についてはチンプンカンプンでした。
Processingを通じて数学を楽しく学習できたことで、より数学、生物に関する好奇心が湧いてきました。
この作品はチューリングパターンに、癌細胞のパターン形成モデルを組み込んで作品にしたものです。
作品のテーマは「バグ」です。システム開発では「バグ」は根絶すべきですが、
ジェネラティブアートは意図的にバグを引き起こし、意図しない美しさを探るという
変わった手法があります。
生命現象における癌はバグのようなものなので、そのモデルを取り入れて芸術の一部にしてしまうというのがアイデアとして面白いと思いました。
Bug
今回はかなり難しく、ボリュームのある内容でしたが、最後まで読んでいただきありがとうございました！
参考文献
Comments