Houdini初心者のためのローレンツ方程式の数値解導出とビジュアライズ

CPSLab Advent Calendar 2019 20日目の記事です。

みなさんは プロシージャル（手続き的） な3Dモデリングの経験があるでしょうか。この記事はプロシージャルモデリングの手法でローレンツ方程式の計算の過程を解説します。

プロシージャルな手法

ググれば腐るほど詳細が出てきますが、ソフトウェアエンジニアの方々ならば「プロシージャルな手法」という言葉に魅了的な響きを感じるものがあるのではないでしょうか？
あまりピンとこなかった方にピンとこないであろう解説をすると、 （数学的な意味も含む）関数の組み合わせ によってモデルを生成するという手法です。

これを実用的なものに落とし込んだものが、プロシージャルな手法でマテリアル（テクスチャ）を作成するSubstance Designer、そして今回紹介・使用するHoudiniです。

Houdiniとは？

プロシージャルモデリングを恐らく最も積極的に行っているソフトウェアであり、水や煙、炎におけるビジュアルエフェクトについて強力なソリューションを有しているソフトウェアです。
試用については無料であり、開発元のSideFXのページからApprentice版を無期限かつ無料利用することができます（制限はありますが、今回は全く問題にならないため省略します）。

インストールすると最上位ライセンスであるHoudini FXなど全部のライセンスバージョンがインストールされますが、起動するのはHoudini Apprenticeです。

今回の目標

次のようなモデルが出力されることをゴールとしています。

よく目にするローレンツ方程式の図ですね。今回はローレンツ方程式が表現する状態を求めることまでを目標にします。ローレンツ方程式についてはググればたくさん出てきますが、難しいことはさておいて以下の方程式で与えられます。

Houdiniことはじめ

まず、Houdiniの細かいUIはさておくとして、プロシージャルモデリングを実施するにあたりHoudiniはノードベースのエディタを利用します。その中で今回利用するノードは大きく分けて SOP(Surcace OPerator 、 Wrangle 、SOP Solver の3種類です。また、Houdiniのプログラミング言語として VEXpression を利用します。
恐ろしいことにHoudiniではこれ以外にもノードの種類があり、プログラミング言語においてはHScript(deprecated)やPython、C++を使う機会も存在しますが、今回は使いません。

SOP

まず先にSOPから説明します。Surface OPeratorの通り、サーフェスに対して行う処理です。が、頂点や辺、点に対してもオペレーションが行われます。これらはそれらを作成する、あるいは変更を加えるものがあります。また、出力や状態を司るものまで多種多様揃っています。まるでRxのようです

基本的に、ルートのノードでは作業を行わず、 Geometry ノード以下で行います。この関係はディレクトリ構成とよく似ています…

ノードエディタ（デフォルトでは画面右下）にフォーカスした状態でTabキーを押下し、geoまで入力するとGeometry SOPが候補に上がります。今回はこれを利用します。
Enterを2回押下し、ノードを配置します。
そしてgeo1ノードをダブルクリックし、ノードの中に入ります。これはサブウィンドウ右上の戻る進むボタンからいつでも上の階層あるいは下の階層に戻ることができ、さながらディレクトリのようです。

たとえば正多面体を生成するPlatonic Solid SOPです。こちらはgeo1ノードに入った状態でふたたびTabキーを押下し、platonicまで押下すると候補が出てくるので、やはり同じように配置します。

このSOPのふるまいはデフォルトの画面右上のアトリビュートエディタから変更可能です。たとえば、デフォルトではSolid TypeがTetrahedron（正三角錐）ですが、画像ではDodecahedron（十二面体）に変更しています。この他にもPlatonicではありませんが、サッカーボールやUtah Teapodなどが選択できます（Utah Teapodは有名ですよね！）

では、Utah TeapodからSOPの強力な機能を見ていきましょう。Platonic Solid SOPにおけるUtah Teapodは（恐らく）NURBSサーフェスの集合ですが、これをポリゴンに変更してみましょう。
手法はおそろしく簡単で、Platonic SolidからConvert SOPノードに接続するだけです。Convert SOPを作成し、Platonic SolidからConvertへ接続し、Convertノードの一番右側、マウスホバー時に青色になる部分をクリックすればポリゴン化されたUtah Teapodが出現するはずです。

あっというまです…このとき、Convert SOPにフォーカスした状態でアトリビュートエディタのU, Vの値を変更すると、それぞれの方向についてポリゴンがどれほどの密度になるのかを変化させることができます。

このとき、ノードを削除あるいはバイパス（ノード左の黄色になる部分を有効）にすれば瞬時にそのオペレーションは無視された状態のオブジェクトが生成されます。

また、数学的な側面も存在します。IsoSurface SOPでは、Implicit Surface（陰関数曲面） 注:正式な和訳を知りません… の描画が行えます。

デフォルトでは球が描画されます。アトリビュートエディタのImplicit Functionを見ると、球の描画に使われるx^2 + y^2 + z^2 - 1 = 0の左辺がそれとなく表現されていますね！このままでは不便なポリゴンなので、Remesh SOPを用いるとメッシュが再構成されます。

ほんの片鱗ですが、これがHoudiniの強力さのあらわれです。ほかにもたくさんの種類のSOPが存在し、それらは非破壊的に行われます…つまり、もとのオブジェクトに変更は一切かかっていないということです（いつでももとに戻せます！）

そして、重要なSOPに、Group SOPが存在します。

これは要するにpointやedgeを選択したものを保存しておくSOPですが、この後のSOPなどでこのGroupを指定してそこだけ効果を適用する、などという器用なことが可能になります。この場合入力したサーフェスのpoint番号1から100までのポイントをGroup1として保存しています。

WrangleノードとVEXPression(VEX)

WrangleノードはPoint Wrangle, Attribute Wrangleなど各種ありますが、どれもVEXpressionというプログラミング言語を利用したSOPです。種類によりますが、たとえばAttribute Wrangleはアトリビュートの調整を行うSOPです。

注：このままでは見づらいので、上記の図では線をポリゴンチューブにするWire SOPを利用しています

このサンプルのために、まずはLine SOPを作成し、アトリビュートのLengthを10に設定します。このノードをResampleノードに接続し、Lengthを0.01に設定します。そしてAttribute Wrangleに接続し、次のコードをアトリビュートエディタ内のコードエディタに記述します…

@P.x = 0;
@P.z = float(@ptnum) / 100;
@P.y = sin(radians(@ptnum));

ソフトウェアエンジニアの方へ…セミコロンでの文の終端から、floatへの型キャスト、関数呼び出しなど、見慣れたオペレータ群が揃っていますね！

セミコロンは文の終端にないとエラーになります。ご注意！

しかし、暗黙的に利用されているアットマーク開始の変数たちにフォーカスする必要がありますね。まず、@Pが何なのかについて解説する必要があるでしょう。これはWrangleが処理しているポイントのひとつを指しており、x, y, zのプロパティを持っています。そして@ptnumは現在のポイントの番号を指しています。

これらをなんとなく理解するために、Attribute Wrangleにフォーカスし、メインの画面からGeometry Spreadsheetを開きます…そこは地獄のようです

落ち着いて見ていきましょう。@Pの持つ値が一覧として表示されています。これらをアトリビュートと呼び、これらの値に代入することでアトリビュートを変更しています。また、アトリビュートは任意にジオメトリに挿入することができます。

f@amount = @ptnum;

fは浮動小数点であるfloatを表しており、float型としてamountを定義し@ptnumの値を代入しています。これをAttribute Wangleに代入すると、Geometry Spreadsheetにamountの項目が追加されていることに気がつくでしょう。これは値を保存し、後々のWrangleで利用したいときに便利です。

SOP Solver DOP

SOP Solver DOPは特殊で、アニメーションにおいて（画面下のタイムラインに関することです）、前のフレームを参照して現在のフレームを作成するDOPとなっています。

…いきなりDOPという言葉が出てきましたね。もう面倒になっているのであまり解説しませんが、Dynamic Operator nodeの略です。文字通りダイナミックに変化するものに対して使えるオペレータですが、私はこれ以外のDOPを知りません…

たとえばsin波をアニメーションしてみましょう。先ほどのネットワークからWrangleを削除し、次のようにSOP Solverを配置します。

そして、SOP Solverにダブルクリックし、SOP Solverの中に入ります。

SOP Solverのネットワークは至極単純です。前のフレームを表すPrev_Frameノード、そしてSOP Solverにつないだ入力4つのノードです。

今回はPrev_Frameを主に利用します。Prev_FrameにAttribute Wrangleを繋ぎ、次のコードを挿入します。

@P.x = 0;
@P.y = sin(radians(@Frame*10 + 500 * @ptnum / 1000));
@P.z = float(@ptnum) / 100;

そして、ウィンドウ左下の再生ボタンをクリックします。するとサイン波が横に移動している（ように見える）でしょう。

青色の部分が計算（キャッシュ）済みの部分です。ここの部分については計算なしで再生できますが、より複雑な計算をSOP Solverで利用するとメモリが枯渇しがちです。

もし計算が重くて、動かなくなってしまった場合、Escキーを押下することで計算をキャンセルできます。覚えておくと便利です。

さあ、ローレンツ方程式を計算しよう！

1. Add SOPで始点を作成

Add SOPは何かと便利なものですが、ここではローレンツ方程式を計算するに当たり、始点をひとつ登録するだけのことをします。
アトリビュートエディタでPointタブを開いた状態で、始点を0.1, 0.1, 0.1にし、チェックボックスを有効化して点を登録します。

1. last Groupを作成する

ここでは何も選択せず、group nameをlastにします。これで後ほどlast groupを利用することができます。今の所、中身は空っぽですが…

1. Solver SOPを作成する

Solver SOPを作成し、中に入りPrev_FrameにAttribute Wrangleを作成します。このとき、WrangleのGroupをlastに設定します。これは"last" Groupにのみこの処理を行うという意味です。この効果は後ほど。

そして、Wrangleは次のようにコードを記述します。

float dx(vector Point) {
    float p = 10.0;
    return p*(Point.y - Point.x);
}

float dy(vector Point) {
    float r = 28.0;
    return Point.x*(-Point.z + r) - Point.y;
}

float dz(vector Point) {
    float b = 8.0/3.0;
    return Point.x * Point.y - b * Point.z;
}

float dt = 0.01;

int newpointnumber = addpoint(geoself(), @P + set(dx(@P), dy(@P), dz(@P)) * dt);
setpointgroup(geoself(), "last", @ptnum, 0);
setpointgroup(geoself(), "last", newpointnumber, 1);

関数の宣言と実装、返り値についてはここで語るまでもないでしょうが、念の為解説をします。

float dx(vector Point) {
    float p = 10.0;
    return p*(Point.y - Point.x);
}

この関数dx()はvector - つまりx,y,z3つの座標(float)を持つ型を受け取り、その値についてfloat型の値を返す関数です。中身の計算については、かつて上に示したローレンツ方程式の一式目とほぼ同様なことが伺えます。
ローレンツ方程式については2式目、3式目についても同様の関数を作成します。 - これがこの処理の中核をなします。

ですが、まだいくつかのことについても解説が必要ですね。
set()関数は与えられた引数の数にもとづいて vector2, vector, vector4, matrix型の値を返す関数です。今回は3つの引数が与えられているため、vector型の値が返ってきます。

geoself()関数は処理している自分自身のジオメトリ番号を示しています。この場合ジオメトリは1つしかありませんから、0です。そのためgeoself()の代わりに0でも差し支えありません。

addpoint()関数はジオメトリに文字通り点を追加し、追加された点の番号を返却します。返却された点番号はint型です。

ここで最も肝要なのは、setpointgroup()関数です。これはpointをどのグループに属させるのかという選択の関数であり、引数は対象となるジオメトリ番号（*注:この関数は現状、0以外を受け取りません)、グループ名、ポイント番号、そして選択されている1かされてないか0をとります。

ここで@ptnumは現在処理しているポイントを指しており、これは"last"グループから非選択とします。そして、新しく生成されたポイントであるnewpointnumberを選択としています。

この"last" Group指定は重要です - このWrangleは最初に"last"グループに対して処理するという宣言をしました。"last" Groupが継続的に変化することによって、最後のポイントに対してのみ処理を行えるのです。

1. 経過確認

ここまでの処理を終えたら、一旦上の階層にまで戻り、SOP Solverを可視状態にし、再生ボタンを押します。
ローレンツ方程式が求まっていく様子が確認できると思います。

tips: アニメーションが短すぎてよくわからない！と思ったら、Solver SOPにフォーカスし、アトリビュートエディタからSub Stepsを10にまで増やしてみましょう。1フレームにつき10回計算が行われるようになり、見かけ上光速になります。また、アニメーションの左下Global Animation OptionsからEndを1000程度にし、Applyすると1000フレームまでアニメートされます。

そして、ある程度まで計算されたローレンツ方程式の結果を見てみましょう。うまくいけば次のようになっているはずです…

点群ですね！これは点群なので、「それらしい見た目」とはなっていません。点と点を結び、線にする必要があります。そこで、ふたたびAdd SOPの出番です！SOP Solverに新たにAdd SOPを接続し、アトリビュートエディタからPolygonsタブ>By Groupを選択すると、線分が作成されます。

線分になったからにはやりたい放題です - ここにPolywire SOPをつないでもよいし、Copy To Point SOPの第2入力に入れて、Sphereなどを第一引数に入れて点をもっと大きく描画することも可能です - これが私達のゴールです。

ですが、もう少しあがいてみましょう。ここからはHoudini特有の「お作法」です。

1. お作法

記述したWrangleの中に直書きの数値があったことに気づきましたか？これらをひとつのCONTROLLERというノードにまとめて一括で管理できるようにします。

AddやSolver SOPと同じ階層にNullノードを作成します。名前をCONTROLLERとします。

アトリビュートエディタのヘッダにある歯車アイコンをクリックし、Edit Parameter Interfaceを選択します。すると次のウィンドウが表示されます。

左側のペインからFloatを選択し、ドラッグ・アンド・ドロップで右側のペインに移動します。これを4回繰り返しましょう。

そして右側ペインから上から順にp, r, b, dtと命名します。このとき、dtは値が大きすぎると数値解が発散してしまうため、Rangeにチェックを入れて0.001から0.01とし、両側の鍵のボタンを押します。

すべてが済んだら"Accept"で変更を保存しましょう。

次のようになっているはずです…

これにそれぞれ、p=10, r=28, b=2.6666666667を挿入します。dtは0.01にします。

そして、SOP SolverのWrangleコードを次のように変更し、コードエディタの右側にあるスライダーと+アイコンを一緒にしたようなボタンを押します。すると、コード下部にパラメータが追加されます。

float dx(vector Point) {
    float p = ch("p");
    return p*(Point.y - Point.x);
}

float dy(vector Point) {
    float r = ch("r");
    return Point.x*(-Point.z + r) - Point.y;
}

float dz(vector Point) {
    float b = ch("b");
    return Point.x * Point.y - b * Point.z;
}

float dt = ch("dt");

int newpointnumber = addpoint(geoself(), @P + set(dx(@P), dy(@P), dz(@P)) * dt);
setpointgroup(geoself(), "last", @ptnum, 0);
setpointgroup(geoself(), "last", newpointnumber, 1);

一旦CONTROLLERノードまで戻り、pの値を右クリックして「Copy Parameter」を実行します。

ふたたびSOP SolverのWrangleノードに戻り、パラメータPの値を一旦殻にしてから右クリックし、今度は「Paste Relative References」をクリックします。

するとch("../../../../CONTROLLER/p")が挿入されます。まさにrelative referencesです。これはCONTROLLER(Null)ノードのパラメータを参照しています。同様に他のパラメータについても同じことを行いますが、面倒なのでこれをコピペし、

ch("../../../../CONTROLLER/r")
ch("../../../../CONTROLLER/b")
ch("../../../../CONTROLLER/dt")

と言った具合に推定できるものはとっとと書いてしまいましょう。さいわいなことに、補完が効いてくれます。

もう一度実行してみましょう…何も変わりません。では、CONTROLLERからパラメータのスライダーを色々いじってみて、ローレンツ方程式がどのように変化するのかを観察してみましょう…これからは遊びの時間です！自由にローレンツ方程式の係数を変えて、様々な形に表情を変えていく様を見ましょう！

補足

説明が完全に入り込む隙間がありませんでしたが、最後の解説としてVEXについて説明せねばならないでしょう。「なぜローレンツ方程式をSOP Solverを経由して解かなければならなかったのか？」
答えはVEXが完全に並行に実行されるプログラミング言語なためです。これはある点あるいはエッジの前後を参照できないということを指しており、これがローレンツ方程式のビジュアライズにとっては非常にネックとなってしまいます…ローレンツ方程式のビジュアライズには、前のポイントの情報が必要なのですから。そのため、今回のローレンツ方程式のビジュアライズはVEXのパワーを活かしきれていないということです…

そう、VEXはとても強力な言語です。意識せず並行処理が行われ、SIMDなどの支援による高速な実行が期待できます。ですから、可能な限りVEXで完結するような処理を考えるべきです。

なお、無視していましたがPythonのノードが用意されており、Pythonでとっととローレンツ方程式を描写することもできます。今回はもうやりませんが。何度も言いますがVEXで書けることはVEXで書くべきなのです。

おわりに

この動画に助けられました: https://www.youtube.com/watch?v=saA6-edb-OE

Houdiniは3Dデザインを主としており、上記の特性を利用してアルゴリズミックなデザインも可能になります…フィボナッチ数を用いて花をそれらしく表示したり、マンデルブロ集合を3Dに拡張したものを表示したり…またその特性から無機物と自然物のモデリングも得意です。そして、忘れないでおいてほしいのが、Houdiniは炎や水、破壊の表現も得意です…これらの組み合わせはとても強力なものばかりです。Houdiniワールドは広い。この荒野をもっと旅する気分になったなら、Appertanceライセンスでいつまでも冒険を続けることが可能です…

ようこそ、プロシージャルモデリングの世界へ！