読者です 読者をやめる 読者になる 読者になる

追実装: Projective Dynamics (SIGGRAPH 2014)

SIGGRAPH 2014 論文 開発 CG

この記事の概要

  • 現在用いられる物理エンジンの多くは position-based dynamics (PBD) [Muller et al. 2007] を基礎としている *1
  • 近年提案された projective dynamics [Bouaziz et al. 2014] は PBD より優れた性質を持っているらしい
  • そこで projective dynamics の追実装を行った上で、その将来性について検討する

f:id:yuki-koyama:20170206150935g:plain
Projective dynamics による布のシミュレーション

なお projective dynamics は以下の記事でも取り上げたことがあります。今回は前回よりさらにマニアックな内容になっています。

PBD については以下の記事で取り上げています。

アルゴリズム解説

準備:Implicit Euler Integration (Section 3.1)

Projective dynamics は implicit Euler integration (後退オイラー積分法) に基づく手法です。したがって数値的に安定な性質があります。

まず運動方程式に implicit Euler integration を適用すると、時刻 t_{n} における状態 *2 と時刻 t_{n+1} における状態の関係式を得ることができます。

f:id:yuki-koyama:20170206153601p:plain

時刻 t_{n} における状態を既知として時刻 t_{n+1} における状態を未知とすると、この関係式は連立方程式とみなすことができます。これを解くのが目標です。

ここで Martin et al. [2011] が提案した式変形を用いることで、この方程式はエネルギー最小化問題 (最適化問題) として定式化することができます。

f:id:yuki-koyama:20170206153848p:plain

この目的関数は、

  • 慣性エネルギーに関する項
  • 弾性エネルギー *3 に関する項

の和となっています。慣性に関する項は単純な quadratic ですが、弾性エネルギーに関する項は通常非線形であり *4、このエネルギー最小化問題を高速に解くのは難しいです。

Non-Linear Elasticity (Section 3.2)

そこでこの論文では、非線形な弾性エネルギーの扱いについて新しい考え方を導入します。

通常の連続体力学における超弾性体 *5 の扱いにおいては、弾性変形状態 \mathbf{q} は歪み (strain) \mathbf{E}(\mathbf{q}) *6 として表現され、これと材料の歪みエネルギー密度関数 (strain energy density function) *7 を用いることで、弾性エネルギー W(\mathbf{q}) が計算されます。

この論文では弾性エネルギー W(\cdot) を、

  • 歪みに関する制約を充たす変形状態がなす多様体 (constraint manifold): \mathbf{E(\cdot) = \mathbf{0}}
  • 二つの変形状態間の距離関数: d(\cdot, \cdot)

という、二つの要素で定義し直します。より具体的には、

\displaystyle{
W(\mathbf{q}) = \min_{\mathbf{p}} \left\{ d(\mathbf{q}, \mathbf{p}) + \delta_{\mathbf{E}} (\mathbf{p}) \right\}
}

と弾性エネルギーを定義し直しています。ここで

\displaystyle{
\delta_{\mathbf{E}}(\mathbf{p}) = \left\{ \begin{array}{ll}
    0 & (\mathbf{E}(\mathbf{p}) = 0) \\
    + \infty & (\text{otherwise})
  \end{array} \right.
}

であり、\mathbf{p}多様体のうちで \mathbf{q} に最も近い変形状態を探すための補助変数です。このように表現すると、弾性エネルギーの計算は

  1. 現在の変形状態を constraint manifold に投影 (projection) する *8
  2. 現在の変形状態と投影後の変形状態の距離を計算する

の二つの処理を行うことで実現できるということになります。

この処理を高速に行うために、projective dynamics では距離関数 d(\cdot, \cdot)\mathbf{p}, \mathbf{q} に対して quadratic に定義することにします。つまり

\displaystyle{
d(\mathbf{p}, \mathbf{q}) = \| \mathbf{A} \mathbf{p} - \mathbf{B}\mathbf{q} \|_{F}^2
}

の形式 *9 で表せるとします。

なお quadratic な距離関数に制限してしまうことのデメリットについては、

The constraint nonlinearity (also known as geometric nonlinearity) is already taken care of by the projection on the constraint set, so the distance metric can be kept simple, trading general material nonlinearity against efficiency and robustness.

と説明されています。つまり、材料の非線形性 (material nonlinearity) は表現できなくなってしまうものの、制約の非線形性 (constraint nonlinearity) は考慮されているから良いだろうということのようです。

Projective Implicit Euler Solver (Section 3.3)

上述の弾性エネルギーの定義を用いると、エネルギー最小化問題として表現された運動方程式

  • Local solve: 投影 (projection) によって \mathbf{p} を探す
  • Global solve: \mathbf{p} を固定した上で \mathbf{q}_{n+1} に関して最小化問題を解く

の二つの処理を繰り返すことで解くことができるようになります。繰り返しの回数を増やすほど真の解に収束していきます。この考え方は global/local alternating optimization または block coordinate descent と呼ばれます *10

Local solve は制約に対して並列化が可能です。具体的に何を計算するかは後で (論文では Section 5) 説明されます。

Global solve は目的関数が \mathbf{q}_{n+1} に対して quadratic であるため、目的関数を \mathbf{q}_{n+1}微分した値が \mathbf{0} になるという性質を用いると、線形方程式として解くことができます。さらにこの線形方程式は、係数行列 *11 が時間不変である *12 ため、最初に一度だけ行列の分解を行っておけば、各タイムステップでは高速に計算することが可能となります。

連続体ベースの制約条件 (Section 5)

この論文では projective dynamics で扱える連続体ベースの制約として

  • Edge strain 制約 (Liu et al. [2013] の手法に類似)
  • Triangle strain 制約 (Chao et al. [2010] や Wang et al. [2010] の手法に類似)
  • Tetrahedron strain 制約
  • 面積保存制約
  • 体積保存制約
  • 例示ベース弾性制約 (Martin et al. [2011] や Koyama et al. [2012] の手法に類似)
  • Bending 制約

などをまとめて提案しています。

これらの連続体ベースの制約は、既存の多くの幾何ベースの制約 *13 と異なり、メッシュの解像度の変化や非一様性に対してロバストに動作するという特徴があります。

その他の制約条件 (Section 6)

この論文では、連続体ベースの制約以外の制約として

  • 位置制約 (ある頂点を引っ張る、空間に固定する、などの用途に使用)
  • 衝突制約

などを説明しています。

さらに、projective dynamics は PBD の一般化であるため、おそらく PBD で用いられる制約 *14 は全て組み込むことが可能です (要確認)

実装

既存のオープンソース実装

EPFL の Computer Graphics and Geometry Laboratory によって開発された実装がオープンソースとして公開されています。

今回実装するにあたって、クラス構造などを中心に ShapeOp による実装を参考にしました。

使用ライブラリ

以下のライブラリを使用しました。

  • Eigen (線形システムソルバー *15、ベクトル・行列クラス)
  • libigl (三次元モデルの読み込み、法線方向の計算)

並列化

C++11 の std::thread を用いてマルチスレッド化を行いました *16。具体的には local solve 部分は制約毎に並列化が可能で、global solve 部分は各座標軸 ( x, y, z) 毎に並列化が可能です。

制約

今回は布のシミュレーションを実装するために最低限な

  • Triangle strain 制約
  • 位置制約

の二つを実装しました。布のシミュレーションとしては

  • Edge strain 制約
  • Bending 制約
  • 衝突制約

なども実装すると良いと思いますが、今回は扱いません。

実行環境

以下の結果は MacBook Pro (Retina, 13-inch, Mid 2014) で実行しています。CPU は 3 GHz Intel Core i7 です。リアルタイム描画ではなく連番画像を出力する実装になっています。

結果

結果1:異なる解像度

この論文で projective dynamics のために提案された連続体ベースの制約群は、メッシュの解像度の変化や非一様性に対してロバストに動作するのが特徴です。以下の動画では異なるメッシュの解像度でも (皺の詳細度の違いこそあれど) 似た結果が得られていることを確認しています。これは従来 *17 の PBD では実現できなかった性質です。

この動画では全ての例で反復回数を 15 回に統一しています。

以下に計算時間を示します。リアルタイム計算と呼べるのが 15 [ms / frame] 程度までと考える *18 と、概ね 5k DoFs まではリアルタイムで動くといって良さそうです。

#triangles #dofs initialize [ms] local [ms / iteration] global [ms / iteration] overall [ms / frame]
64 123 0.096 0.147 0.111 4.011
256 435 0.263 0.191 0.138 5.108
1,024 1,635 0.949 0.306 0.180 7.553
4,096 6,339 4.245 0.704 0.374 16.78
16,384 24,963 28.63 2.309 1.589 60.28
65,536 99,075 256.7 9.619 9.344 291.8

結果2:異なる反復回数

Projective dynamics の最大の特徴は、反復回数を増やせば増やすほど正しい *19 解に収束していくことです。以下の動画では反復回数をだんだんと増やしたときに挙動が収束していくことを確認しています。これも従来 *20 の PBD では実現できなかった性質です。

この動画では全ての例でメッシュ解像度を #triangles = 16,384 に統一しています。計算時間はほぼ反復回数に比例しています。この例では 10 回程度反復すれば十分な印象です。解像度を高くしたり、硬い材質を表現するために制約の重みを大きくしたりすると、より多くの反復が必要になると思います。

Projective Dynamics は PBD に取って代わるものになるか?

Projective dynamics の論文を読む限り PBD よりもずっと良さそうに見えますが、将来どちらが市民権を獲得することになるかは今のところ判断できないというのが結論です。理由は以下の 4 点です。

理由1:2014 年以降に発表された PBD の拡張手法が結構良さそう

Projective dynamics の論文が投稿された 2014 年 1 月の段階では、PBD には

  • 連続体ベースの制約が扱えないため、結果がメッシュに依存してしまう
  • 反復回数に依存して結果が変わってしまう

という欠点がありました。今回の論文ではこれらの欠点を解消することに成功しています。

しかしその後、PBD で連続体ベースの制約が扱えるようになった [Muller et al. 2014] (要確認) り、反復回数に依存して結果が変わってしまう問題を解消する (要確認) 手法 XPBD が提案されたり [Macklin et al. 2016] してきています。これらによって projective dynamics を採用する利点は小さくなっているように思えます。

理由2:Projective dynamics は PBD に比べ計算時間がかかる

実は projective dynamics の local solve は PBD の反復計算とほぼ同等の処理を行っています。したがってこの部分の計算コストはほぼ同等と考えられます。

一方で projective dynamics には global solve の処理が存在します。この処理によって反復回数に依存して結果が変わってしまう問題を解消できる反面、PBD よりも計算コストが上がってしまっています。あくまで目安ですが、今回の実装ではおよそ 1.5 倍から 2 倍程度、global solve によって全体の計算コストが上昇したと考えられます。

Local solve は GPU を駆使することで格段に速くなる余地がありますが、global solve はこれ以上の並列化は困難です。したがって、どうしても projective dynamics は PBD よりも大幅に計算時間がかかってしまうことになります。

追記 (2017/02/09):Global solve を GPU で並列化する論文が Wang [2015] によって提案されていました。ご指摘ありがとうございます。この論文では global solve にあたる線形方程式の計算を GPU 上で Jacobi 法を使って解く方法を説明しており、さらに Chebyshev 準反復法と呼ばれる手法を組み合わせることで projective dynamics の大幅な高速化を行っています。

  • Huamin Wang. 2015. A chebyshev semi-iterative approach for accelerating projective and position-based dynamics. ACM Trans. Graph. 34, 6, Article 246 (October 2015), 9 pages. DOI: http://dx.doi.org/10.1145/2816795.2818063

理由3:Projective dynamics は PBD に比べ柔軟性に欠ける

Projective dynamics では global solve の係数行列を事前に分解しておくことで実行時には高速に方程式を解くことができます。しかし、制約が変わる場合はこの計算コストの高い行列分解を再計算する必要が生じてしまうことから、PBD に比べ柔軟性に欠けるという問題があります。

例えばパーティクルベースの流体シミュレーションでは、粒子間の制約が (ほぼ) 毎ステップ、動的に変わっていきます。柔軟性に欠ける projective dynamics では、毎ステップ行列分解をするのを避けるために、matrix-free CG 法を用いて枠組みに変更を加えた上で流体シミュレーションを実現する必要があります [Weiler et al. 2016]。一方で PBD では行列分解を必要としないため、ストレートフォワードに流体シミュレーションを実現することができます [Macklin and Muller 2013]。

理由4:Projective dynamics も 2014 年以降だんだんと進化してきている

PBD だけでなく、projective dynamics もまた 2014 年以降少しずつ拡張手法が提案されてきています。様々な材質の超弾性体 *21 をシミュレーションすることができるようになったり [Narain et al. 2016; Liu et al. 2016]、ハード制約 *22 が扱えるようになったり [Narain et al. 2016] してきています。

(一応) 結論

Projective dynamics は計算コストを小さく保ったままより物理的に正確な表現を実現する方向で発展していることもあり、映像制作にはかなり向いている印象です。一方で PBD はその計算スピード・柔軟性・近年の拡張手法の存在から、ゲームや VR などのインタラクティブコンテンツでは依然として選ばれ続ける気がします。

最後に

本記事に誤りや分かり難い点などを見つけた際には、些細な点でも構いませんので、ご連絡いただけると嬉しいです。

References

*1:Bouaziz et al. [2014] によれば PhysX、Havok Cloth、Maya nCloth、Bullet は PBD に基づいているそうです。その他にも Maya nHair や Houdini も PBD だと思います。

*2:ここでは m 個の頂点が存在している世界を考え、それら頂点の位置 \mathbf{q} \in \mathbb{R}^3 と速度 \mathbf{v} \in \mathbb{R}^3 によって状態が記述されると考えます。

*3:より正確には内力のポテンシャルエネルギーですが、この論文では特に弾性エネルギーにフォーカスを絞って議論を進めていきます。

*4:線形な弾性エネルギーを定義することも可能ですが、コンピュータグラフィクスで求められる大変形などを表現するには不適切だということだと思います。

*5:超弾性体とは(大まかには)理想的な弾性特性を持つ物質のことで、微小変形ではなく大変形の解析に用いられます。

*6:グリーン歪み (Green strain) がよく使われます [Martin et al. 2011]。

*7:Saint Venant-Kirchhoff model や Neo-Hookean model など材料のモデル毎に定義されます。

*8:これが projective dynamics の名前の由来だと思われます。

*9:なお  \mathbf{A} = \mathbf{B} = \mathbf{I} とすると単純なユークリッド距離を考えていることになります。

*10:近年 Liu et al. [2013] によって使用されました。

*11:ここでは線形方程式 \mathbf{Ax} = \mathbf{b}\mathbf{A} の部分を指します。

*12:これが偶然なのか必然なのかは謎です。

*13:従来の PBD における制約は幾何的なものがほとんどで、メッシュの解像度などが変化すると挙動が大きく変わってしまうという問題がありました。

*14:Shape matching 制約などが含まれます。

*15:今回は疎なエルミート行列を分解するために Eigen::SimplicialLDLT を使用しました。

*16:感覚的にはマルチスレッド化によって 2 倍程度処理が速くなった印象です。

*17:この論文が投稿された時点より前を指します。

*18:もちろん、実際にゲームや VR コンテンツに組み込むにはこれよりもずっと小さな計算コストに抑える必要があります。

*19:物理的に正確という意味ではなく、エネルギー最小化問題の解として正解に近いという意味です。

*20:この論文が投稿された時点より前を指します。

*21:例えば Saint Venant-Kirchhoff model や Neo-Hookean model も projective dynamics の枠組みでシミュレーションできるようになったそうです。

*22:通常の PBD や projective dynamics では制約は全てソフト制約として扱われています。