HoudiniのVEXを用いた雪の結晶の生成シミュレーションの一例 | An example of snowflake growth simulation using VEX in Houdini

こちらの記事は Houdini Advent Calendar 2021 14日目の記事です。

はじめに

　去年に引き続き今年もHoudini Advent Calendarに参加しました！
　都内のCGスタジオでshotを作ったり、ツールを設計開発したりするTalkyと申します。

　本記事の内容は以前参加した案件に使う素材の試作品を手直した物となります。発注元から「雪の結晶を降らしたい」という要望があったので、でも単にIllustratorで雪の結晶っぽいイラストを描いてHoudiniやMayaで立体化させるのは面白くないから、遊び半分で本物のシミュレーションの論文を参考にしました。
　論文に基づいたシミュレーションであるため、少々堅苦しい解説が多くなりますが、参考になれば幸いです。

作業環境

• Windows 10 Pro 64bit
• Houdini FX 19.0.383 (コアの部分はAttribute Wrangleしか使わないので、H17.5以降は動作確認済みです)

1. 資料検索

　Googleで「Houdini Snowflake Paper」を検索すると優先順位の上に来る結果はEntagmaさんの「Creating A Procedural Snowflake」と堀川淳一郎先生の「Three-Dimensional Snowflake」二つの動画になります。
　2つのチュートリアルはそれぞれC. A. Reiter, A local cellular model for snow crystal growth (2005)とJ. Gravner, D. Griffeath, Modeling snow crystal growth: a three-dimensional mesoscopic approach (2009)の理論部分を参考として使われました。
　Houdini界隈で代表的な方のチュートリアルはどちらも素晴らしい内容ですが、シミュレーション用のパラメータの数と実行速度は私が欲しい物と少し違うので、自力で資料となる論文を探すことに決めました。
　そこで、チュートリアルではなく元となる論文を色々調べた結果、J. Gravner, D. Griffeath, Modeling Snow Crystal Growth II: A mesoscopic lattice map with plausible dynamics (2008)はパラメータの数が丁度いいので、こちらの論文¹の理論部分を参考することにしました。

2. 生成モデル

　もし§1.で言及した二つのチュートリアルを試作したことがある方でいれば、Houdini上の作業はさほど変わることがないので、元の論文もしくはこの章で私の解説をベースにして自力で作成してもいかがでしょう。

2.1. 生成フィールド

　雪の結晶の生成フィールドは辺の長さが統一された正六角柱が平面上に満遍なく配列させると前提条件とします。それぞれの正六角柱が一つの格子と考え、幾つの分量[mass]と状態を持ちます。
　単位時間毎に中心となる格子と近傍の6つの格子の間に各分量が予め与えられた環境パラメータを使って計算処理し、格子の状態は各分量の数値と環境パラメータに従って切り替えます。これらの過程を繰り返して格子の結晶状態をシミュレーションすることにします。

2.2. 分量と状態

　全ての格子に「結晶分量」[Crystal Mass]、「拡散分量」[Diffusive Mass]、「境界分量」[Boundary Mass]と3つの分量を持ちます。
　そして格子は「雪の結晶の一部」（結晶格子）、「結晶周りの蒸気」（蒸気格子）、「結晶と蒸気の境界」（境界格子）といずれの状態になっています。
　Houdini上でこちらの要件を再現する場合は格子がPoint、分量がPoint Attribute、状態がPoint Groupを使います。
　なお、こちらの分量と状態の名前は英語で表記した物以外、私は記述の便宜上に付けた名前です。

2.3. 環境パラメータ

　環境パラメータは各格子の間に分量の変動や特定の状態になるかどうかの閾値になる。言い換えると、雪の結晶を作り上げるためのルールになります。パラメータを変えることによって雪の結晶の生成シミュレーション結果も変わります。
　各パラメータの説明は下表と§2.4を合わせてご覧ください。

パラメータ	説明	推奨範囲
Rho(ρ)	蒸気格子と境界格子が持つ拡散分量の初期値	[0.3, 0.9]
Beta(β)	一部の境界格子が結晶に切り替えるための境界分量の閾値	[1.05, 3.0]
Alpha(α)	一部の境界格子が結晶に切り替えるための境界分量の閾値、θと一緒にknife-edgeを構成する	[0.0, 1.0] θより大きい
Theta(θ)	一部の境界格子が結晶に切り替えるための拡散分量の閾値、αと一緒にknife-edgeを構成する	[0.0, 1.0] αより小さい
Kappa(κ)	境界格子が凍結時に持つ拡散分量が結晶分量へ変わる凍結係数、残りは境界分量へ変わる	[0.0, 0.2]
Mu(μ)	境界格子が融解時に持つ境界分量が拡散分量へ変わる気化係数	[0.0, 0.2]
Gamma(γ)	境界格子が融解時に持つ結晶分量が拡散分量へ変わる昇華係数	[0.0, 0.01] μより遥かに小さい
Sigma(σ)	拡散分量に影響を与えるノイズ	[0.0, 0.00005]

2.4. 結晶の過程

　結晶の過程は単位時間毎に全て格子に対して「拡散」[Diffuse]、「凍結」[Freezing]、「付着」[Attachment]、「融解」[Melting]の順番で計算処理を行います。
　ノイズ パラメータSigma(σ)が0ではない時に「融解」後に「拡散ノイズ」[Noise]の過程を追加します。
　全ての過程の処理が終わりましたら次の単位時間のループに入ります。こうした単位時間のループを1万回から10万回を繰り返します。Frame Rangeを10万にするのは現実ではないのでSolver SOPのSub StepsでFrame毎のループ回数を増やします。

2.4.1. 拡散[Diffuse]

　拡散の過程は「境界格子」と「蒸気格子」に対して処理を行います。
　1) 対象格子と近傍の6つの格子が持つ拡散分量の総和を求めます。もし近傍の格子が結晶格子である場合、結晶格子の拡散分量は0.0のため、対象格子自身の拡散分量に置き換えます。
　2) 拡散分量の総和を格子総数の7で割って対象格子の拡散分量に再代入します。

2.4.2 凍結[Freezing]

　凍結の過程は「境界格子」に対して処理を行います。
　1) 対象格子が持つ拡散分量に凍結係数Kappa(κ)を掛けて、格子が持つ結晶分量に累加します。
　2) 残りの拡散分量は境界分量に累加します。（係数は1-κになります）
　3) 対象格子の拡散分量を0.0にします。

2.4.3 付着[Attachment]

　付着の過程は「境界格子」に対して処理を行います。この過程はこのシミュレーションの肝心な部分であります。
　対象格子の近傍格子の状態によって3つの分岐条件があります。それぞれの分岐で「付着成功」かどうかを判定して最後は各分量の計算処理を行います。

2.4.3.1 近傍格子の結晶格子個数は1と2の場合

　対象格子が持つ境界分量が閾値Beta(β)より大きいか又は等しい場合、「付着成功」になります。

2.4.3.2 近傍格子の結晶格子個数は3の場合

　1) 対象格子が持つ境界分量が1.0より大きいか又は等しい場合、「付着成功」になります。
　2) もしくは、近傍格子が持つ拡散分量の総和が閾値Theta(θ)より小さい、且つ対象格子が持つ境界分量が閾値Alpha(α)より大きいか又は等しい場合、「付着成功」になります。

2.4.3.3 近傍格子の結晶格子個数は4以上の場合

　対象格子の分量状態と関係なく、「付着成功」になります。

---

上記の3つの分岐で「付着成功」と判定された対象格子が持つ境界分量が結晶分量に累加し、境界分量を0.0にします。格子の状態を「結晶格子」に切り替えます。

2.4.4 融解[Melting]

　融解の過程は「境界格子」に対して処理を行います。付着の過程で「付着成功」にならない境界格子はこの過程を経て各分量が少し変更されて次の単位時間のループに入ります。
　1) 対象格子が持つ境界分量に気化係数Mu(μ)を掛けて、拡散分量に累加します。残りの境界分量は再代入します。（係数は1-μになります）
　2) 対象格子が持つ境界分量に昇華係数Gamma(γ)を掛けて、拡散分量に累加します。残りの結晶分量は再代入します。（係数は1-γになります）

2.4.5 拡散ノイズ[Noise]

　拡散ノイズの過程はノイズ パラメータSigma(σ)が0ではない時に「蒸気格子」と「境界格子」に対して処理を行います。
　対象格子が持つ拡散分量にSigma(σ)を掛けて、累加もしくは累減します。それぞれの確率は0.5になります。

---

　以上は雪の結晶の生成シミュレーションのモデルになります。ここからはHoudiniで作業をします。

3. Houdiniでの作業

3.1 Controller Nodeの作成

　最初にフィールドやシミュレーションに必要なパラメータを設定できるController Nodeを作ります。
　空のGeometry Objectの中にNull SOPを作成し、Edit Parameter Interfaceで下表にあるパラメータを追加します。

パラメータ	変数名	タイプ	範囲（Range）	初期値	説明
Grid Row	grid_row	Integer	[50, 1000]	300	生成フィールドの層数（大きさ）
Grid Interval	grid_interval	Float	[0.1, 10.0]	0.5	Pointの間隔
Core Row	core_row	Integer	[0, 10]	0	最初の結晶核の層数（大きさ）
Rho (ρ)	rho	Logarithmic Float	[0.0, 1.0]	0.8	環境パラメータ（拡散分量の初期値）
Beta (β)	beta	Logarithmic Float	[1.0, 3.0]	1.9	環境パラメータ（付着の閾値）
Alpha (α)	alpha	Logarithmic Float	[0.0, 1.0]	0.004	環境パラメータ（付着の閾値）
Theta (θ)	theta	Logarithmic Float	[0.0, 1.0]	0.001	環境パラメータ（付着の閾値）
Kappa (κ)	kappa	Logarithmic Float	[0.0, 0.2]	0.05	環境パラメータ（凍結係数）
Mu (μ)	mu	Logarithmic Float	[0.0, 0.2]	0.015	環境パラメータ（気化係数）
Gamma (γ)	gamma	Logarithmic Float	[0.0, 0.01]	0.0001	環境パラメータ（昇華係数）
Sigma (σ)	sigma	Logarithmic Float	[0.0, 0.00005]	0.0	環境パラメータ（拡散分量のノイズ）
Random Seed	random_seed	Float	[0, 10]	任意	乱数生成用パラメータ
Sub Steps	sub_steps	Integer	[10, 100]	50	Solver SOPのSub Steps制御用

　なお、私は各分量を可視化するために3つのColor Rampも作成しましたが、こちらのパラメータは必ずしも作成する必要がないです。

3.2 フィールドの生成

　Controller NodeのGrid Rowパラメータは初期値の300のままだと途中でViewportがPointに埋め尽くされるので、ここで一旦10に設定します。

　Add SOPでフィールドの中心となるCenter Pointを一つ作成します。

---

　次はAttribute Wrangle SOP(Run Over: Detail)を1つ作成し、周りのPointを生成します。Add SOPの出力とAttribute Wrangle SOPの第一入力を接続します。
　フィールドの形について、今回のシミュレーションは正六角形の雪の結晶なので、四角のフィールドを作ると角辺りのPointは生成されても最後までほぼ使われない場合が多いです。余計な計算量を削減するため、最初から正六角形のフィールドを生成します。最長辺が同じの正方形と正六角形の面積で比較すると約35%の計算量を節約できます。
　更に、こちらのコードを使って生成されたPointの順番は中心から外側まで昇順になるので、Geometry Spreadsheetでシミュレーションの数値を確認する時も分かりやすいです。未検証の仮説ですが、Point番号の走査も中心から始まるので、この記事の事例ではVEXの実行速度にも良い影響があるかもしれません。

　生成コードはforループを使ったシンプルな物です：

generate_hexagonal_points(Detail)

// Get center point position
vector center_pt_pos = point(0, "P", 0);

// Get parameters
float interval = chf("grid_interval");
int row_total = chi("grid_row");
int column_total = 6;

// Generate hexagonal points
for (int row = 1; row < row_total; row++) {
    for (int column = 0; column < column_total; column++) {
        vector pt_pos, next_pt_pos;
        pt_pos.x = sin(2.0 * $PI * column / column_total) * interval * row + center_pt_pos.x;
        pt_pos.y = center_pt_pos.y;
        pt_pos.z = cos(2.0 * $PI * column / column_total) * interval * row + center_pt_pos.z;

        int next_column = column + 1;
        next_pt_pos.x = sin(2.0 * $PI * next_column / column_total) * interval * row + center_pt_pos.x;
        next_pt_pos.y = center_pt_pos.y;
        next_pt_pos.z = cos(2.0 * $PI * next_column / column_total) * interval * row + center_pt_pos.z;

        for (int inter = 0; inter < row; inter++) {
            vector inter_pt_pos;
            inter_pt_pos = (next_pt_pos - pt_pos) * (float(inter) / row) + pt_pos;
            addpoint(0, inter_pt_pos);
        }
    }
}

　最初はAdd SOPで作成されたCenter Pointの座標を取得します。Pointは1つしかないので、Point番号は0になります。
　必要な外部パラメータはController NodeからNode Parameter経由で取得します。

　次はFor文の2重ループを使って階層数row・頂点番号column・格子の間隔intervalと三角関数の組み合わせで六角形の頂点の座標pt_posと次の頂点の座標next_pt_posを計算します。
　注意が必要のは、Center Pointはすでにあったので、rowの初期値は1からとなります。

　計算で得たpt_posとnext_pt_posを使って、もう一つのFor文で頂点を含めてのPointを生成します。頂点と次の頂点の間にあるPointの数はちょうどrowと一致するので、rowを上限としてそのまま使います。

3.3 Pointの初期化

　Attribute Wrangle SOP(Run Over: Points)で全てのPointに対して初期化処理を行います。第一入力は前のgenerate_hexagonal_pointsAttribute Wrangle SOPの出力と接続します。

---

　Pointの初期化は主に3つの部分に分かれます：

• 各Point（格子）のNear Points（近傍格子）を取得し、Point Attributeとして各Pointに格納
• 各PointのPoint AttributeとPoint Groupの初期値を設定
• 最初の結晶核となるPointの範囲を取得し、Point AttributeとPoint Groupの初期値を設定

　Point Attributeは4つあります：

Point Attribute	データ型	説明
boundary_mass	float	境界分量
crystal_mass	float	結晶分量
diffusive_mass	float	拡散分量
near_points	int[]	各格子自身のPoint番号と近傍格子のPoint番号

　Point Groupは3つあります：

Point Group	説明
boundary	境界格子
crystal	結晶格子
vapor	蒸気格子

---

　VEXコードも簡単な物です：

initialize_points(Points)

// Get parameters
float interval = chf("grid_interval");
float rho = chf("rho");
int core_row = chi("core_row");

// Get near points
int near_points[] = nearpoints(0, @P, interval * 1.01);

// Set attribute and group to points which have 6 near points
if (len(near_points) == 7) {
    i@group_vapor = 1;
    i[]@near_points = near_points;
}

// Set attribute and group to all points
i@group_boundary = 0;
f@boundary_mass = 0.0;
f@diffusive_mass = rho;

// Get core points max count
int core_pt_max = (core_row * (core_row + 1) / 2) * 6;

// Set attribute and group to core points
for (int core_pt = 0; core_pt < core_pt_max + 1; core_pt++) {
    setpointgroup(0, "crystal", core_pt, 1);
    setpointgroup(0, "vapor", core_pt, 0);
    setpointattrib(0, "crystal_mass", core_pt, 1.0);
    setpointattrib(0, "diffusive_mass", core_pt, 0.0);
}

　まずは必要な外部パラメータはController NodeからNode Parameter経由で取得します。
　rhoは拡散分量の初期値を示すパラメータになります。

　// Get near pointsの部分は全てのPointの座標に対して、Point間隔であるintervalを半径として近傍のPoint番号を取得します。小数精度の問題で稀に取得できない問題回避するため間隔intervalに1.01を掛けます。
　取得したPoint番号はnear_points配列に代入します。

　ここで注意するのはnearpoints関数は座標と半径でPoint番号を取得するため、近傍Point以外に対象Point自身の番号も取得されます。
　従って、フィールド一番外側を除く、シミュレーションに使えるPointのnear_points配列の長さは7となります。
　こちらのPointにnear_points配列をPoint Attributeとして格納し、蒸気格子を示すPoint Group vaporを設定します。

　// Set attribute and group to all pointsの部分は全てのPointに対し共通のPoint AttributeとPoint Groupを設定します。
　フィールド一番外側のPointはシミュレーションのループ処理は回しませんが、近傍Pointとして拡散の処理に参加するので、拡散分量を与えます。
　境界格子を示すPoint Groupは単位時間のループ毎に変わるのでSolver SOP内で詳細な処理を書きますが、ここで予め空のPoint Groupのみ作成します。
　境界分量boundary_massも初期値の0.0を与えるます。
　拡散分量diffusive_massの初期値rhoを与えます。

　// Get core points max countの部分は最初の結晶核となるPoint番号の最大値を取得します。
　フィールドは中心から生成したのおかげで、結晶核の層数core_rowとPoint最大値の関係は下記の単純な数式で計算できます。VEXはsummation functionがないため、普通に展開して計算します。

max = 6\left( \sum_{k=0}^n k \right)

　for文を使って全ての結晶核Pointcore_ptに対してPoint AttributeとPoint Groupを設定します。Point番号の最大値までループが必要ので、条件式の上限にご注意ください。
　結晶格子を示すPoint Group crystalに設定し、同時に蒸気格子Group vaporから取り除きます。
　結晶分量crystal_massの初期値1.0を与えて、拡散分量diffusive_massを0.0に設定します。

　こちらのVEXに関して、入れ子を使えばもっと綺麗に効率よく書けることはもちろん存じておりますが、説明しやすいため、あえて平にしました。
　なお、この記事のVEXコードは実行速度重視のため、@attributeと@groupを多用してます。分かりにくくて申し訳ございません。

---

　Node Infoで全てのPoint AttributeとPoint Groupが正しく作成されたことを確認します。

　Pointの初期化が終わったので、次はシミュレーションの部分に入ります。

3.4 Solver SOPでシミュレーション

　Solver SOPを1つ作成します。忘れない内にController NodeからSub Stepsのパラメータを取得します。

　Solver SOPの内部に入って、Attribute Wrangle SOP(Run Over: Points)を1つ作成します。第一入力のみPrev FrameNodeの出力と接続します。
　Groupをcrystalに設定して、結晶格子Pointに対して処理を行います。
　この過程は結晶格子の周りにある蒸気格子を境界格子に切り替えます。

set_boundary_group(Points)

foreach (int near_point; i[]@near_points) {
    if (inpointgroup(0, "vapor", near_point) == 1) {
        setpointgroup(0, "vapor", near_point, 0);
        setpointgroup(0, "boundary", near_point, 1);
    }
}

　まずはPoint Attribute near_pointsに格納されたPoint番号を取り出します。
　全ての近傍格子に対して、もしvaporのPoint Groupに入ってあれば、vaporを解除し、boundaryに切り替えます。

3.4.1 拡散の処理

　Attribute Wrangle SOP(Run Over: Points)を1つ作成します。Groupをboundaryとvaporに設定して、境界格子と蒸気格子Pointに対して拡散の処理を行います。

diffusion(Points)

float diffusive_mass_cur = 0.0;

foreach (int near_point; i[]@near_points) {
    float diffusive_mass_prev;

    if (inpointgroup(0, "crystal", near_point) == 1) {
        diffusive_mass_prev = f@diffusive_mass;
    }else{
        diffusive_mass_prev = point(0, "diffusive_mass", near_point);
    }

    diffusive_mass_cur += diffusive_mass_prev;
}

diffusive_mass_cur /= 7.0;

f@diffusive_mass = diffusive_mass_cur;

　拡散分量の変数名はcurとprevのsuffixで現在Frameの状態と前Frameの状態区別します。
　Point Attribute near_pointsに格納されたPoint番号を取り出します。
　§2.4.1で説明した通り、もし近傍格子が結晶格子の場合、対象格子自身の拡散分量を取得します。結晶格子以外はそれぞれの格子が持つ拡散分量を取得します。
　全ての拡散分量を累加して、最後は7で割って、対象格子に再代入します。

3.4.2 凍結の処理

　Attribute Wrangle SOP(Run Over: Points)を1つ作成します。Groupをboundaryに設定して、境界格子Pointに対して凍結の処理を行います。
　Controller Nodeから凍結係数Kappa(κ)のパラメータをNode Parameter経由で取得します。

freezing(Points)

float kappa = chf("kappa");

f@boundary_mass += (1.0 - kappa) * f@diffusive_mass;
f@crystal_mass += kappa * f@diffusive_mass;
f@diffusive_mass = 0.0;

　§2.4.2で説明した通り、各境界格子に対して拡散分量に凍結係数kappaを掛けて、格子が持つ結晶分量に累加します。
　もう一度拡散分量に(1.0 - kappa)を掛けて、境界分量に累加します。
　最後は拡散分量を0.0に設定します。

3.4.3 付着の処理

　Attribute Wrangle SOP(Run Over: Points)を1つ作成します。Groupをboundaryに設定して、境界格子Pointに対して凍結の処理を行います。
　Controller Nodeから閾値Beta(β)、Alpha(α)、Theta(θ)のパラメータをNode Parameter経由で取得します。

attachment(Points)

// Get parameters
float beta = chf("beta");
float alpha = chf("alpha");
float theta = chf("theta");

// Count crystal in near_points and calculate total diffusive_mass
int near_crystal_count = 0;
float diffusive_mass_total = 0.0;
foreach (int near_point; i[]@near_points) {
    if (inpointgroup(0, "crystal", near_point) == 1) {
        near_crystal_count += 1;
    }
    diffusive_mass_total += point(0, "diffusive_mass", near_point);
}

// Test if boundary can attachment
if (((near_crystal_count == 1 || near_crystal_count == 2) && f@boundary_mass >= beta)
    || (near_crystal_count == 3 && (f@boundary_mass > 1 || (diffusive_mass_total < theta && f@boundary_mass >= alpha)))
    || (near_crystal_count >= 4)){

    f@crystal_mass += f@boundary_mass;
    f@boundary_mass = 0.0;

    i@group_boundary = 0;
    i@group_crystal = 1;
}

　// Count crystal in near_points and calculate total diffusive_massの部分はnear_pointsからPoint番号を取得します。
　近傍格子の中から結晶格子の個数を数えて変数near_crystal_countに累加します。
　対象格子と近傍格子が持つ拡散分量の総和を変数diffusive_mass_totalに格納します。

　// Test if boundary can attachmentは近傍結晶格子の個数と拡散分量の総和を§2.4.3で説明した条件で、
　「付着成功」となると対象格子が持つ境界分量を結晶分量に累加し、境界分量を0.0に設定します。
　格子の状態もboundaryからcrystalに切り替えます。

3.4.4 融解の処理

　Attribute Wrangle SOP(Run Over: Points)を1つ作成します。Groupをboundaryに設定して、境界格子Pointに対して融解の処理を行います。
　Controller Nodeから気化係数Mu(μ)と昇華係数Gamma(γ)のパラメータをNode Parameter経由で取得します。

melting(Points)

float mu = chf("mu");
float gamma = chf("gamma");

f@diffusive_mass += mu * f@boundary_mass + gamma * f@crystal_mass;
f@boundary_mass = (1.0 - mu) * f@boundary_mass;
f@crystal_mass = (1.0 - gamma) * f@crystal_mass;

　§2.4.4で説明した通り、境界分量に気化係数muを掛けて、結晶分量に昇華係数gammaを掛けて、拡散分量に累加します。
　境界分量に1.0 - muを掛けて、境界分量に再代入します。結晶分量に1.0 - gammaを掛けて、結晶分量に再代入します。

3.4.5 拡散ノイズの処理

　Attribute Wrangle SOP(Run Over: Points)を1つ作成します。Groupをboundaryとvaporに設定して、境界格子と蒸気格子Pointに対して拡散の処理を行います。
　Controller NodeからノイズパラメータSigma(σ)と乱数パラメータRandom SeedをNode Parameter経由で取得します。

noise(Points)

float sigma = chf("sigma");
float random_noise = rand(54.52 * @ptnum + 14.26 * @SimFrame + 89.64 * chf("random_seed") + 83.16);

if (random_noise < 0.5) {
    f@diffusive_mass += sigma * f@diffusive_mass;
}else{
    f@diffusive_mass -= sigma * f@diffusive_mass;
}

　§2.4.5で説明した通り、乱数によって格子が持つ拡散分量にsigmaを掛けて、0.5の確率で拡散分量に累加または累減します。
　乱数を取得する時に注意が必要のはSolver SOPのSub Stepsが設定されたので、Sub Steps毎に違う乱数を取得したい場合は通常のフレーム数を表すグローバル変数@Frame($F)ではなくSolver SOPなどDOP専用の変数@SimFrame($SF)を使うべきです。

　ノイズの過程はノイズパラメータSigma(σ)が0の時は省略できますので、Switch-If SOPを1つ作成します。ExpressionのEnableに[Sigmaパラメータのパス] != 0.0を記入します。私の場合はスクリーンショットの通りch("../../../../CONTROLLER/sigma") != 0.0になります。
　なお、Switch-If SOPが存在しないHoudiniのバージョンではSwitch SOPのSelect Inputに同じ内容を記入すれば大丈夫です。

　融解過程のmelting Attribute Wrangle SOPからの出力を左の入力（False）に接続して、ノイズ過程のnoise Attribute Wrangle SOPからの出力を右の入力（True）に接続します。

---

最後にOutput SOPを1つ作成し、結果を出力します。シミュレーションの主要な過程はこれで終わりになります。

3.5 シミュレーション結果の処理

3.5.1 シミュレーション結果の確認

　Geometry Objectに戻って、Controller NodeのGrid Rowパラメータを初期値の300に戻します。各パラメータの数値（Random Seedは任意で大丈夫です）はこのようになります。

　Solver SOPを選択してGeometry Spreadsheetでデータを確認します。
　§3.1で提示した各パラメータの初期値と同じように設定された場合は数値は上図のようになります。違った場合は上流の過程に設定ミスがあると思いますので、各Attribute Wrangle SOPのInput・Group・Run Over・VEX・Parametersが正しく設定されるかどうかを確認してください。

3.5.2 シミュレーション結果をキャッシュ

　初歩の確認が終わりましたらFile Cache SOPとNull SOPを1つずつ作成して、シミュレーションの結果を20フレームほどキャッシュを取ります。

　ViewportのカメラをTopに切り替えて、G又はHキーでViewportのカメラをフィールドの中心に移動します。
　ご覧の通り、シミュレーション結果はPoint AttributeとPoint Groupとして保存されたので、Pointの見た目は何も変わらないです。

　結果を直感で確認するため、簡易的な可視化Nodeを作成します。

3.5.3 シミュレーション結果の可視化

　Attribute Wrangle SOP(Run Over: Points)を1つ作成して、Ramp(Color) Parameterを3つ作成します。

visualization(Points)

if (i@group_crystal == 1) {
    // Crystal point
    v@Cd = chramp("crystal_color", clamp(f@crystal_mass / 2.0, 0.0, 1.0));
}else if (i@group_boundary == 1) {
    // Boundary point
    v@Cd = chramp("boundary_color", clamp(f@boundary_mass / 2.0, 0.0, 1.0));
}else if (i@group_vapor == 1){
    // Vapor point
    v@Cd = chramp("vapor_color", clamp(f@diffusive_mass, 0.0, 1.0));
}else{
    // Unused point
    v@Cd = {1, 0, 0};
}

　格子状態のPoint Groupによってそれぞれの分量のPoint Attributeを取得し、PointのCd色を設定します。
　crystal_massとboundary_massの値は大体[0.0, 2.0]にあるので、2.0で割ります。
　clamp関数を使わないくても大丈夫ですが、Ramp ParameterのPosition値域である[0.0, 1.0]を越えた場合、見た目に影響があるので、値の上下限を設定します。

　Viewportに戻りまして、綺麗な雪の結晶が見えてきましたね。同じ状態のPoint Groupに入っても分量のPoint Attribute値の違いも一目で確認できます。

　もっと成長した雪の結晶を見たい場合はキャッシュを120フレームまで取ることをおすすめします。結果は下図のようになります。

4. 雪の結晶を立体化

　ここからはシミュレーションの結果を更に処理を加えて立体のGeometryを作成します。

4.1. 結晶を作成

　Blast SOPを1つ作成します。crystalPoint Group以外のPointとPoint Groupを削除します。

　Tube SOPを1つ作成し、Columnsを6設定します。これを結晶格子の単位構造として使えます。
　Radius Scaleに[grid_intervalパラメータのパス] / sqrt(3)を記入します。私の場合はch("../CONTROLLER/grid_interval") / sqrt(3)になります。
　数学的な計算の過程を省けますが、六角柱の半径は近傍Pointとの間隔の$1/\sqrt{3}$倍になります。

　Copy to Points SOPで単位構造を全てcrystal Pointに複製します。

　ViewportのカメラをPrepsに戻して、雪の結晶が立体化されたことを確認します。
　このままだと雪の結晶がただの板になるので、crystal_massの数値を使って雪の結晶の表面に凹凸を作ります。

　Attribute Wrangle SOP(Run Over: Points)を1つ作成して、Blast SOPとCopt to Points SOPの間に入れます。

crystal_mass_to_scale(Points)

if (f@crystal_mass < chf("threshold")) {
    removepoint(0, @ptnum);
}else{
    float scale_y = pow(f@crystal_mass * chf("scale"), chf("pow"));
    v@scale = set(1.0, scale_y, 1.0);
}

　threshold、scale、powのNode Parameterを追加します。
　thresholdはcrystal_massの下限閾値になります。普段は0のままで大丈夫ですが、稀に生成されたcrystal_massが小さい薄い結晶を削除する時に使います。
　scaleとpowでY軸方向のスケールとスケールの強弱をつけるために使われます。デフォルト値は0.7と2に設定します。

　計算されたscale_yをset関数でfloat3に変換してscale Point Attributeに代入します。
　Y軸方向のスケールのみを制御したいので、均等スケールのpscaleではないことをご注意ください。

　これで雪の結晶の表面に凹凸を作られました。§3.5.3で可視化された明度が高い箇所の単位構造が長くなり、逆に明度が低い箇所の単位構造が短くなります。

　最後はAttribute Delete SOPとGroup Delete SOPを1つずつ作成して、Pとscale Point Attribute以外のPoint AttributeとPoint Groupを全て削除します。
　データ量を軽減のため、私は普段に使われてないAttributeとGroupの即時に解放する癖がありますが、こちらのDelete Nodeを追加しなくても下流の作業に影響がないです。

4.2. 雪の結晶の最終仕上げ

　使い道によって、このままの形でもいいですが、滑らかなポリゴンが欲しい場合はよく使われるPolygonをVDBに変換して形を整えて、Polygonに戻すこともできます。
　こちらの処理の流れは人それぞれですが、ここは最も簡単なVDB from Polygons SOP、VDB Smooth SDF SOP、Convent VDB SOPの流れでスムーズします。

　キャッシュと取る時や他のツールへ書き出し時のデータ量を節約するため、PolyReduce SOPとPack SOPを1つずつ作成します。

4.3. キャッシュと書き出し

　最後はFile Cache SOPとNull SOPを1つずつ作成して、形を整えた雪の結晶のGeometryのキャッシュを取って、流れの終点を示すOUTPUTのNull SOPで終わります。

　後は必要に応じてAlembicにしたりHDAにしたりするのは皆さんにお任せします。

5. おわりに

　雪の結晶のシミュレーション、いかがでしょうか。前人の知恵である学術論文のモデルを拝借して、パラメータ制御によって本物っぽい雪の結晶を生成してみました。
　解説の文章を書くと長くなりますが、実際にHoudini上の作業はとてもシンプルなNodeと簡単なVEXの組み合わせだけです。是非皆さんも作ってみて、色んなパラメータの組み合わせを試してください。

　そう言えば、Advent Calendarの参加を決めた時にTwitterで「去年に作ったものだけど、少しだけVEXを書き直した。」と呟いたな、あれは噓だ。1年前に自分が書いたVEXコードを振り返ったら、色々耐えられなくて全てのVEXを書き直しました。
　VEXの最適化と自分の作業環境がH17.5からH19.0にアップグレードしたことのおかげで、キャッシュを取る時間は以前作業した時の1/10に短縮することもできました。

　何か面白い題材があったら、来年もHoudini Advent Calendarに参加したいと思います！
　最後までご覧いただき、誠にありがとうございました！

　Talky Ren　2021/12/14

Reference

• J. Gravner, D. Griffeath, Modeling Snow Crystal Growth II: A mesoscopic lattice map with plausible dynamics, Physica D: Nonlinear Phenomena 237 (3), 385-404, 2008.

---

1. 参考文献は上記となります。 ↩