動機

Houdiniには剛体をアニメーションする機能として、BulletとかHoudini RBDとかODEという剛体エンジンを使ったDOPノードが用意されているのですが、所詮ブラックボックスなので、中身を自由に改造して機能を拡張することは難しかったりします。

そこで、VEXというHoudiniに内蔵されている言語を使って、剛体のソルバを自作してみました。この記事では慣性テンソルを使った正確な剛体の回転をVEXで書く方法について解説したいと思います。

GitHubに*.hipファイルをアップしていたので適宜参考にして下さい。
https://github.com/nobuyuki83/houdini_vex_rigidbody

剛体の運動方程式

剛体とは回転と並進の自由度だけで記述される（つまり弾性変形しないような）質点の集合です。ここでは、基準配置からの回転がと並進がそれぞれ、直行行列$R\in SO(3)$と、３次元ベクトルで$P$を用いて表現されるとしましょう。つまり、基準配置で$p_0$の位置にあった点は、剛体変形後に$p=R p_0+P$の位置にあるとしましょう。

剛体の角速度については２通りの異なる定義があるのですが、ここでは以下のように角速度ベクトル$\Omega$を定義します:

\dot{R} = R\tilde{\Omega}

$\Omega$の上にチルダが乗っていますが、これはベクトル$\Omega$から作られる反対称(anti-symetric)な行列を意味しています。要するには外積を行列表記にしたものですね。詳しくはWikipediaの記事を参考にしてください:Wikipedia: 交代行列

ここでは詳しくは解説しませんが、この角速度ベクトルは次のように時間発展します（外力＝０とします）。

\dot{\Omega} = I^{-1} (\tilde{\Omega} I \Omega)

ここで出てくる$I$は3x3の行列で、慣性テンソルと呼ばれるものです。HoudiniのRBD ObjectではRotational Stiffnessという名前のパラメータです。慣性テンソルを使うと回転のエネルギー$W$は

W = \frac{1}{2}\Omega^T I \Omega

のように書くことができます。例えば、ある軸$u$を考えたときに、$u^T I u$が大きければ、「その軸回りに回転させるのには大きなエネルギーが必要」なので「軸$u$の周りには回転しにくい」ということです。

慣性テンソルは必ず正定値対称なので、互いに直行する慣性主軸とう軸を持ちます。慣性主軸と$\Omega$が並行だと、$\Omega$と$I\Omega$が並行になるので、外力が無いと$\dot{\Omega}=0$となり、角速度は変化しません。一定の速度で一定の向きに回転します。逆に、慣性主軸と$\Omega$が並行でないと、角速度が時間変化します。これを歳差運動(Procession)といいます。

これは私の個人的な意見ですが、この歳差運動が有ると、剛体の回転の「動きのリッチさ」が増します。なので、慣性テンソルを考えて正確にアニメーションすることは結構重要なんじゃないかと思っています。

VEX プログラミング

アトリビュートの生成

XYZ軸に直行するような、X辺、Y辺、Z辺の長さがそれぞれ、１，２，３である直方体がランダムな初期角速度で回転することを考えましょう。この直方体の慣性テンソル$I$は次のように書けます。

I = \left[\begin{array}{lll}
13 & 0 & 0\\
0 & 10 & 0\\
0 & 0 & 5\\
\end{array}\right]

次のようにVEXを書くと、頂点のアトリビュートとして、慣性テンソルや質量、角速度、初期姿勢を生成できます。

v@W = set(1.0*rand(@ptnum*2)-0.5, 
          4.0*rand(@ptnum*4)-2.0, 
          4.0*rand(@ptnum*8)-2.0); // initial Angular velocity
p@orient = {0,0,0,1}; // initial orientation in quaternion
vector ex = set(1,0,0);
vector ey = set(0,1,0);
vector ez = set(0,0,1);
matrix3 Ix = outerproduct(ex,ex);
matrix3 Iy = outerproduct(ey,ey);
matrix3 Iz = outerproduct(ez,ez);
3@I = 13*Ix + 10*Iy + 5*Iz; // inertia tensor

角速度の速度

前章の式をそのままVEXに書き起こすと、角速度の速度は以下のように計算できます。

vector dw = invert(3@I)*cross(v@W, 3@I*v@W);

前進オイラー法

前進オイラー法を使うと、姿勢と角速度は次のように更新されます。

@orient = qmultiply(@orient, quaternion(@TimeInc,v@W) );
v@W = v@W + dw*@TimeInc;

おわりに

前進オイラー法は時間ステップが大きくなると不安定になります。GitHubのコードでは４次のRunge-Kutta法を使って時間積分をしていますので、ぜひ参考にして下さい。

参考

理論と実践で学ぶHoudini -SOP&VEX編-
https://www.amazon.co.jp/dp/4862463592/