こちらは、 Houdini Apprentice Advent Calendar 2025 の7日目の記事になります。

Houdiniで生命のパターンを自動生成してみる！ -反応拡散系-

それはある日、ウチの子のスグちゃん (スタンデングヒルヤモリ) をお世話していたとき、「このあまりにも美しすぎる背中の模様はなんなのか…」と疑問に思い、調べてみました。

[ スグちゃんの美しすぎるお姿 ]

なんでも、一般的に自然界にはこのような模様が現れることがしばしばあり、これを「チューリングパターン」というそうです。
言わずと知れたアラン・チューリングがこの原理をはじめて数学的に記述し、20世紀末よりコンピュータの計算能力の急速な発達によってこの分野の研究の発展がもたらされました。

ということで、今回はこのチューリングパターンを実際にHoudini上で実装してみようと思います。ご興味ある方はHIPも添付しますのでそちらも見ていってください！

実装途中で気づいたんですが、Houdini21.0から「reactiondiffusion_block」が実装されており、1から実装しなくてもデフォルトノードで試せるようになってました。SideFXさんお目が高い...
ということで、ここ以降はご興味ある方は覗いてみてください！

1.数式を読み解く

まずは、チューリングパターンを表わす数式から読み解いていきましょう。今回では、比較的簡単な Gray-Scott モデルを使用します。

\frac{\partial u}{\partial t}
  = D_u \nabla^{2} u -uv^2 + f(1-u),

  \frac{\partial v}{\partial t}
  =  D_v \nabla^{2} v + uv^2 - (f + k)v

ここで、f (feed) は供給を表し、k (kill) は消滅を表します。
この数式を一つずつ解説していきます。
まず、左辺のこの部分ですが、

\frac{\partial u}{\partial t}, \frac{\partial v}{\partial t}

下がdtとなっているので、uとvの成分をそれぞれ時間方向に微分するということです。今回はCGソフトでの実装になるので、「1フレームごとにこの計算をやればいいんだな」と理解してしまって大丈夫です。

次に、この部分です。

D_u \nabla^{2} u , D_v \nabla^{2} v 

▽^2と書いてあるのはラプラシアンというもので、関数の勾配を表わすみたいです。ここで『拡散』を表現します。ただ、これは実装としては思ったより簡単なので、身構える必要はないです。

最後に、この部分です。

-uv^2 + f(1-u),

uv^2 - (f + k)v

ここで『反応』を表現します。二つのパラメータ f, k によって u, v の反応をコントロールします。
この式によって、u が活性化すると v が抑制するようになり、f と k のパラメータが良いバランスだと、u と v が偏在的に共存するようになり、結果として縞模様やまだら模様などのパターンが生み出されます。これが、「チューリング不安定性」です。

2.実装

本実装に入る前に、OpenCLに慣れてない方はこちらに目を通すことをオススメします。

まず実装に入る前に、@KERNEL と @WRITEBACK の説明をしていきます。Houdini OpenCL では、コードを複数回実行するためにこの機能を実装しています。

#bind layer &src float
#bind layer &dst float

@KERNEL
{
    @dst.set(@src);
}

@WRITEBACK
{
    @src.set(@dst);
}

@KERNEL のコードで @dst に計算した値を入れ、@WRITEBACK で @src に @dst の値を戻すことで、並列計算時にすでに計算して書き込み済みの値を呼び出すことを回避しています。

---

さて、それでは準備ができたので、数式をもとに Houdini OpenCLCOP で実装していこうと思います。
まず拡散項です。

#bind layer &src_a float
#bind layer &dst_a float
#bind layer &src_b float
#bind layer &dst_b float

@KERNEL
{
    int2 ixy = @ixy;
    
    // 4方向をサンプリング
    float a_c = @src_a.bufferIndex(ixy);
    float a_w = @src_a.bufferIndex((int2)(ixy.x - 1, ixy.y));
    float a_e = @src_a.bufferIndex((int2)(ixy.x + 1, ixy.y));
    float a_n = @src_a.bufferIndex((int2)(ixy.x, ixy.y + 1));
    float a_s = @src_a.bufferIndex((int2)(ixy.x, ixy.y - 1));
    
    float b_c = @src_b.bufferIndex(ixy);
    float b_w = @src_b.bufferIndex((int2)(ixy.x - 1, ixy.y));
    float b_e = @src_b.bufferIndex((int2)(ixy.x + 1, ixy.y));
    float b_n = @src_b.bufferIndex((int2)(ixy.x, ixy.y + 1));
    float b_s = @src_b.bufferIndex((int2)(ixy.x, ixy.y - 1));
    
    // ラプラシアン近似
    float lap_a = (a_w + a_e + a_n + a_s) - (a_c * 4.0f);     
    float lap_b = (b_w + b_e + b_n + b_s) - (b_c * 4.0f);  

    @dst_a.set(lap_a);
    @dst_b.set(lap_b);
}

@WRITEBACK
{
    @src_a.set(@dst_a);
    @src_b.set(@dst_b);
}

a と b、2種類の物質を用意し、それぞれ4方向のサンプリングからラプラシアンを近似します。ちなみにこの状態では a と b どちらも干渉しないですが、アニメーションさせて iteration 数を増やしてみると面白いパターンができあがります。

[ radius=0.1 の円で iteration=49 の結果 ]

次に、反応項を追加し、コードを完成させます。

#bind layer &src_a float
#bind layer &dst_a float
#bind layer &src_b float
#bind layer &dst_b float

#bind parm da float
#bind parm db float
#bind parm feed float
#bind parm kill float
#bind parm dt float

@KERNEL
{
    int2 ixy = @ixy;
    
    // 8方向をサンプリング
    float a_c = @src_a.bufferIndex(ixy);
    float a_w = @src_a.bufferIndex((int2)(ixy.x - 1, ixy.y));
    float a_e = @src_a.bufferIndex((int2)(ixy.x + 1, ixy.y));
    float a_n = @src_a.bufferIndex((int2)(ixy.x, ixy.y + 1));
    float a_s = @src_a.bufferIndex((int2)(ixy.x, ixy.y - 1));
    
    float b_c = @src_b.bufferIndex(ixy);
    float b_w = @src_b.bufferIndex((int2)(ixy.x - 1, ixy.y));
    float b_e = @src_b.bufferIndex((int2)(ixy.x + 1, ixy.y));
    float b_n = @src_b.bufferIndex((int2)(ixy.x, ixy.y + 1));
    float b_s = @src_b.bufferIndex((int2)(ixy.x, ixy.y - 1));
    
    // ウェイトを付けて合成
    float lap_a = (a_w + a_e + a_n + a_s) - (a_c * 4.0f);     
    float lap_b = (b_w + b_e + b_n + b_s) - (b_c * 4.0f); 

    // 反応項
    float reaction = a_c * b_c * b_c;
    
    // Gray-Scottの反応拡散方程式の更新計算
    float delta_a = (@da * lap_a) - reaction + (@feed * (1.0 - a_c));
    float delta_b = (@db * lap_b) + reaction - ((@feed + @kill) * b_c);
    float next_a = a_c + (delta_a * @dt);
    float next_b = b_c + (delta_b * @dt);
    next_a = clamp(next_a, 0.0f, 1.0f);
    next_b = clamp(next_b, 0.0f, 1.0f);
    
    @dst_a.set(next_a);
    @dst_b.set(next_b);
}

@WRITEBACK
{
    @src_a.set(@dst_a);
    @src_b.set(@dst_b);
}

パラメータはこちらを参考にしてください。適切なパラメータでないとパターンが現れないことがあります。

src_a には val=1.0 の白背景、src_b にはお好きな図形やノイズなどを試してみると良いです。
この状態でアニメーションさせてみると、パターンが作られていきます。

以下、他のパラメータを入力した例を挙げます。

[F = 0.03, k = 0.055]

[F = 0.02, k = 0.05]

[F = 0.018, K = 0.045]

付属のHIPには、pigheadの背中にチューリングパターンを出させる作例を載せたので見てみてください。(HIPなどでご質問などある方はDMにてご質問ください)

まとめ

この反応拡散系は、自然界に起きている現象を非常に簡単な方程式で説明できており、自然の神秘に驚かされますね。アラン・チューリング先生が魅了されるのも納得です。
Houdinistとしては、もちろん動物の表面のテクスチャ作成に有用だとは思いますが、その他の用法も考えられそうです。reactiondiffusion_blockを上手に活用していきたいですね！
ということで、本記事は以上になります。よいお年を！

参考

https://kaityo256.github.io/sevendayshpc/day5/index.html
https://www.mls.sci.hiroshima-u.ac.jp/ryo/mathsciA/lecture/Note6.pdf