こちらの記事は Houdini Advent Calendar 2023 19日目の記事です。

はじめに

　2年ぶりに Houdini Advent Calendar に参加しました！
　都内の CG スタジオでエフェクトを作ったり、ツールを設計開発したり、偶にディレクターも務める Talky と申します。

　この記事は xyzdist(), primuv(), prim_normal() などの VEX 関数を使って Scatter SOP で生成された「空間中」の Point に対して Attribute Interpolate SOP で限定的に補間させる方法の一例です。

ご注意

　「限定的」という修飾語を使う理由は、Primitive 単体の最小単位である三角ポリゴンは直交座標系内で同じ平面上にあり、一部ジオメトリ全体に対する変換操作の情報を保持できないからです。
　例えば、X-Z 座標平面上にある三角ポリゴンは Y 軸方向のスケール変形をかけても形は変わらなく、外部情報を経由しない限り変形の結果は保持されないです。（任意の三角ポリゴンも同じく、法線方向のスケールをかけても形は変わりません。この点について文章最後の補足説明をご覧ください。）
　この記事に書かれた方法も Primitive 単体が保持した幾何学的情報に依存しているので、同様に移動と回転以外の変換は動作できますが期待通りの結果は得られません。

動作環境

Windows 10 Pro 64bit
Houdini FX 20.0.506 (Python 3.10)

1. Attribute Interpolate SOP とは

　公式ドキュメンテーションを読んだ方が色々詳細なパラメータが書かれていますが、正直に言うとあのページを読んでノードの動作を理解できる人はかなり少ないと思います。
　よく使われる場面の一つとして、Scatter SOP でジオメトリの表面に Point をばら撒いた時、元のジオメトリが動いたら、Point とジオメトリの相対位置がズレる現象の対策です。言葉だけでは伝わりにくいので、こちらの動画で確認した方が早いと思います。

　実際のノードネットワークはこのようにセットアップされます：

　Time Shift SOP でジオメトリを特定のフレームに固定し、そこから Scatter SOP の Output Attributes タブの Prim Num Attribute と Prim UVW Attribute を ON にして sourceprim と sourceprimuv 2つのアトリビュートを Point に作成します。
　Attribute Interpolate SOP の左入力は生成された動かない Point で、右入力は元の動いているジオメトリと繋ぎます。そうすれば Attribute Interpolate SOP は Prim 番号と Parametric UV 情報† をベースにして、Point を正しい位置に補間（移動）してくれます。

　上記の例以外の使い方もいくつかありますが、調べたい方はすあまさんのブログ記事がおすすめです。

よく Scatter SOP の Relax Iterations オプションが悪さしたとされますが、実は Relax Iterations オプションを OFF にした状態でも、Point はズレること場合もあります。 Scatter SOP の内部アルゴリズムの問題かもしれませんが、確実に対策したい時は Attribute Interpolate SOP を使いましょう。

2. 空間中の Point の場合

　Crag を使うと結果を確認しにくいので、ここからはシンプルな Box をサンプルとして使います。

　ジオメトリ表面上に生成された Point は上記の方法で解決できますが、ジオメトリから Volume (VDB) に変換し、Scatter SOP で空間中に Point を生成する場合、設定を同じにしても、sourceprim と sourceprimuv アトリビュートの値は正しく作成されません。（もちろん Attribute Interpolate SOP で補間を行った場合は想定外の結果になります。）

　そこで、今回の記事のテーマである「空間中」の Point に対して補間する方法を試行錯誤しながら構築します。

2.1. Point に一番近い Prim の情報を生成

　表面上の Point と同様に sourceprim と sourceprimuv 2つのアトリビュートを利用する方針にしたいので、 xyzdist() 関数で Point から元のジオメトリ上に一番近い Point （以降「最短 Point」と呼びます）を探して、hit した Prim 番号と Parametric UV 情報を sourceprim と sourceprimuv に保存します。

　Attribute Wrangle SOP(set_attributes) の1番目の入力に Point のデータ、2番目の入力にフレームが固定された元のジオメトリと繋ぎます。

set_attributes (Points)

// Run Over: Points

int source_prim;
vector source_prim_uv;

xyzdist(1, v@P, source_prim, source_prim_uv);

i@sourceprim = source_prim;
v@sourceprimuv = source_prim_uv;

2.2. 確認するための可視化

　この記事では補間の結果を視認しやすいため、近くの Point 同士の間を繋いだ Polyline を追加し、更に Point と最短 Point に繋ぐ Polyline も別色で表示するようにします。
　実際のセットアップ時に使わない箇所なので、ノードネットワークを確認する時にご注意ください。

　記事のテーマから離れた内容であるため、VEX コードを折りたたみましたが、気になる方はどうぞご確認ください：

可視化 Polyline を追加する VEX コード

generate_polyline (Points)

// Run Over: Points

// Connect near points
float _radius = chf("radius");
int _maxpoints = chi("max_points");

int near_points[] = pcfind(0, "P", v@P, _radius, _maxpoints);
removevalue(near_points, i@ptnum);  // remove ptnum itself from search result

foreach (int near_pt_num; near_points) {
    int near_line = addprim(0, "polyline", i@ptnum, near_pt_num);
    setprimattrib(0, "Cd", near_line, {1.0, 1.0, 1.0});
}

// Add a polyline to hit position
vector hit_pos = primuv(1, "P", i@sourceprim, v@sourceprimuv);
int hit_pt = addpoint(0, hit_pos);

setpointattrib(0, "sourceprim", hit_pt, i@sourceprim);
setpointattrib(0, "sourceprimuv", hit_pt, v@sourceprimuv);

int hit_line = addprim(0, "polyline", i@ptnum, hit_pt);
setprimattrib(0, "Cd", hit_line, {0.5, 0.2, 0.2});

2.3. 自作の sourceprim と sourceprimuv のみで補間した結果

　この2つのアトリビュートのみを使った場合、全ての Point が最短 Point の位置に補間されました。
　もちろん正しい結果にはまだ遠いですが、少なくとも補間した結果はちゃんと元のジオメトリの上に収めることができました。

　ここで足りない情報を考えると、先ず思い付くのは Point から最短 Point までの距離です。
　§2.1 の手順に xyzdist() 関数でアトリビュートの情報を生成する時に最短距離も返り値として出力されるので、その値を Point のアトリビュートとして保持して、補間の後に復元すれば、Point が空間中のどの場所にあったのは再現できると思います。

　但し、xyzdist() 関数の返り値は距離の長さ（float）だけになるので、位置情報として、さらに最短 Point から Point までの方向（vector）が必要です。
　そこで、出番になるのは primuv() 関数です。ジオメトリの Prim 番号と Parametric UV 情報を入力すれば、指定したアトリビュートの補間した値を出力される関数なので、xyzdist() で取得したデータから最短 Point の座標を計算してくれます。

2.4. Point から最短 Point までの位置情報の生成と復元

　§2.1 に作成された Attribute Wrangle SOP(set_attributes) の VEX に少し変更を加えて、最短 Point から Point の位置情報を生成します。

set_attributes (Points)

// Run Over: Points

int source_prim;
vector source_prim_uv;

xyzdist(1, v@P, source_prim, source_prim_uv);

i@sourceprim = source_prim;
v@sourceprimuv = source_prim_uv;

+ vector position_on_prim = primuv(1, "P", source_prim, source_prim_uv);
+ vector source_offset = v@P - position_on_prim;

+ v@source_offset = source_offset;

（ご注意：差分表示のため、先頭に + を記述しましたが、VEX コードをコピー&ペーストする時に取り除いてください）

　この source_offset アトリビュートを利用し、Attribute Interpolate SOP の補間結果を最短 Point から Point の位置に復元します。
　Attribute Interpolate SOP の後に Attribute Wrangle SOP(restore_position) を追加し、中身の VEX を記述します。

restore_position (Points)

// Run Over: Points

v@P += v@source_offset;

2.5. 復元した結果を確認

　Attribute Wrangle SOP(restore_position) の実行結果を確認すれば、Point が「正しい」位置に復元されたことになるのは分かります。
　元のジオメトリを動かしても、ちゃんと空間中の位置が保持されましたね。

　が、

　察しの良い方は気付いたと思いますが、この方法だと、元のジオメトリの Prim が回転した場合 source_offset にその回転情報が反映されず、間違った結果になります。
　試しに元のジオメトリに回転を加えたら、復元結果が崩れたことを確認できます。

　では、Prim の回転情報をどこから取得すればいいと考えると、たぶん法線の方向が一番取得しやす情報で、回転前の Prim の法線方向と回転後の Prim の法線方向から回転行列まは四元数を算出すれば、Prim の回転も復元できます。
　今度は prim_normal() 関数を使って Prim 番号と Parametric UV 情報から Prim の法線方向を取得したいと思います。

2.6. 回転情報の生成と復元

　再び Attribute Wrangle SOP(set_attributes) に戻り、今度は Prim の法線方向を source_prim_normal として Point に保持させます。

set_attributes (Points)

// Run Over: Points

int source_prim;
vector source_prim_uv;

xyzdist(1, v@P, source_prim, source_prim_uv);

i@sourceprim = source_prim;
v@sourceprimuv = source_prim_uv;

vector position_on_prim = primuv(1, "P", source_prim, source_prim_uv);
vector source_offset = v@P - position_on_prim;

v@source_offset = source_offset;

+ v@source_prim_normal = prim_normal(1, source_prim, source_prim_uv);

　Attribute Wrangle SOP(restore_position) の2番目の入力には Attribute Interpolate SOP の右入力と同様に動いている元のジオメトリを入力します。
　下記の VEX で前後法線方向の回転値を計算し、source_offset に回転を加えます。

restore_position (Points)

// Run Over: Points

+ vector cur_prim_normal = prim_normal(1, i@sourceprim, v@sourceprimuv);

+ matrix3 rot_mat = dihedral(v@source_prim_normal, cur_prim_normal);
+ vector adjusted_offset = v@source_offset * rot_mat;

- v@P += v@source_offset;
+ v@P += adjusted_offset;

　私個人は回転行列に慣れたのでこちらを使っていますが、四元数の方が見やすいと思う方は下記のコードを使っても同じ結果になります。（時間計算量的に若干違いますが）

Using Quaternion

vector cur_prim_normal = prim_normal(1, i@sourceprim, v@sourceprimuv);

vector4 rot_quat = dihedral(v@source_prim_normal, cur_prim_normal);
vector adjusted_offset = qrotate(rot_quat, v@source_direction);

2.7. 回転を配慮した復元結果の確認

　今度は元のジオメトリが回転しても、補間後の復元がきちんと行われて、正しい結果になります。

　元のジオメトリを Crag のアニメーションに切り替えても、形状が補間されたことを確認できます。

3. 移動変換と回転変換以外の場合

　記事の最初に書いていますが、今回の Primitive 単体の幾何学的情報に依存する方法は移動と回転しか対応しません。
　参考として、破綻した補間結果も載せます：

　スケールの場合、source_offset に適切なスケール情報を与えることができませんので、最短 Point との距離がずっと同じのままです。
　拡大の場合はスカスカになるが、縮小の場合は反対側からはみ出します。

　せん断（shear）の場合でも、Prim の回転は適用されたものの、やはりスケールは正しく反映できず、一部の場所の歪みは強くなります。

おわりに

　使い道がかなり限られていますが、Attribute Interpolate SOP を使って、空間中の Point の位置を限定的に補間させる方法を1つ提示させて頂きました。

　以前某プロジェクトに参加した時、この問題に直面したけど、色んなところで検索をかけてもなかなか解決方法が見つからなかったため、何時間にかけて自力でこの中途半端な答えにたどり着いた覚えがありました。
　記事を作成時に内容を振り返ったらかなりシンプルな VEX になりますが、せっかくなので当時の試行錯誤を含めて書きました。誰かの力になれば幸いです。

　元々はこの方法を使った作例を1つ作成したいのですが、内容がかなり離れたため、Apprentice の方にもう1つの記事を書く予定になります。

　Talky Ren　2023/12/19

Suppl. 補足説明

Suppl.1. 法線方向スケールの数学的な解釈

　直交座標系の X-Z 座標平面上の三角ポリゴンの各頂点を $(x_0, 0, z_0)$, $(x_1, 0, z_1)$, $(x_2, 0, z_2)$ とします。

　Y軸（法線）方向にスケール $s_y$ を適用する変換行列 $M_{\text{s}}$ があるとして、

{M}_{\text{S}} = \begin{bmatrix}
1 & 0 & 0 \\
0 & s_y & 0 \\
0 & 0 & 1 \\
\end{bmatrix}

　任意の頂点座標に対して変換を適用しても、計算結果を見ればわかると思いますが、変換行列の情報は残されませんでした。

\begin{bmatrix}
1 & 0 & 0 \\
0 & s_y & 0 \\
0 & 0 & 1 \\
\end{bmatrix}
\begin{bmatrix}
x_n \\
0 \\
z_n \\
\end{bmatrix} =
\begin{bmatrix}
x_n \\
s_y \cdot 0 \\
z_n \\
\end{bmatrix} =
\begin{bmatrix}
x_n \\
0 \\
z_n \\
\end{bmatrix}

　X-Z 座標平面ではなく、任意の三角ポリゴンに対してもローカル座標系に変換すれば、同じ結果になります。

Suppl.2. Parametric UV とは

　テクスチャ作業などに使う uv アトリビュートとは別に、各 Primitives が持つ固有な空間の補間情報になります。
　もっと詳しい説明は公式ドキュメンテーションをご覧ください