2017-11-02

Directional Field (方向場) の種類を整理する

論文 CG

Directional field (方向場)*1 や vector field (ベクトル場) についてのメモです。次の論文を参考にしています。

Directional Field Synthesis, Design, and Processing
Amir Vaxman, Marcel Campen, Olga Diamanti, Daniele Panozzo, David Bommes, Klaus Hildebrandt, Mirela Ben-Chen
Eurographics 2016 - State of the Art Reports
https://doi.org/10.1111/cgf.12864

この論文では三次元空間中の surface 上の directional fields (や vector fields) を主に想定して解説しています。

Directional field にはその他にも様々な種類があり、それぞれ様々な名前で呼ばれているようです。また、

Unfortunately, the available terminology in the literature suffers from many inconsistencies [...]

とあるように、時として一貫性のない言葉遣いもされているらしく、注意が必要です。

種類と特徴

代表的なもので次図のようなものが存在します。

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

円は対象の空間中のある点、矢印や線分はその点が持つ情報を表しています。

なお "direction" は方向のみを持つ情報、"vector" は方向と大きさを持つ情報として言葉を使い分けています。また rotationally-symmetric directional field を縮めて RoSy field と呼ぶそうです。前図のうち line field と cross field はそれぞれ RoSy field の一種です。

Vector Field (ベクトル場)

Vector field では空間中の各点が持つ情報は 1 つの vector によって定義されます。例えば scalar field (スカラー場) の勾配を計算することで得られるのも vector field です。

Line Field (線場)

Line field では空間中の各点が持つ情報は 180 度の回転対称性のある 2 つの directions によって定義されます。2-RoSy field とも呼ばれます。Vector field に比べて、magnitude がないだけでなく、前と後ろのどちらを向いているかの区別がありません。各点で定義される値が対称な階数 2 の tensor で表されることから symmetric tensor fields (対称テンソル場) と呼ばれることもあるようです*2。

Cross Field (交差場)

Cross field では空間中の各点が持つ情報は 90 度の回転対称性の関係にある 4 つの directions によって定義されます。4-RoSy field とも呼ばれます。Non-photorealistic rendering (NPR) でも使われます*3。

Frame Field

Frame field では空間中の各点が持つ情報は、必ずしも直交しない 2 つのベクトルと、その 180 度の回転対称性の関係にあるもう 2 つのベクトルによって定義されます。Panozzo et al.*4 によって導入されて以来よく使われているようです。

*1:この訳が一般的かどうかは不明。

*2:Kai Xu, Lintao Zheng, Zihao Yan, Guohang Yan, Eugene Zhang, Matthias Niessner, Oliver Deussen, Daniel Cohen-Or, and Hui Huang. 2017. Autonomous reconstruction of unknown indoor scenes guided by time-varying tensor fields. ACM Trans. Graph. 36, 6, pp.202:1--202:15. URL: http://vcc.szu.edu.cn/research/2017/tfnav/

*3:Aaron Hertzmann and Denis Zorin. 2000. Illustrating smooth surfaces. In Proceedings of the 27th annual conference on Computer graphics and interactive techniques (SIGGRAPH '00). ACM Press/Addison-Wesley Publishing Co., New York, NY, USA, 517-526. DOI=http://dx.doi.org/10.1145/344779.345074

*4:Daniele Panozzo, Enrico Puppo, Marco Tarini, and Olga Sorkine-Hornung. 2014. Frame fields: anisotropic and non-orthogonal cross fields. ACM Trans. Graph. 33, 4, Article 134 (July 2014), 11 pages. DOI: https://doi.org/10.1145/2601097.2601179

2017-03-21

論文読み: Differential Coordinates for Interactive Mesh Editing (SMI 2004)

2004 CG SMI 論文

以下の論文に関するメモです。

Differential Coordinates for Interactive Mesh Editing
Yaron Lipman, Olga Sorkine, Daniel Cohen-Or, David Levin, Christian Roessl, Hans-Peter Seidel
International Conference on Shape Modeling and Applications (SMI) 2004

論文について

SMI 2004 という IEEE 系の国際会議で発表された論文です。この会議は現在では消滅してしまったか改名してしまったかで、そのままは残っていなさそう*1です。

Google Scholar によると 2017 年 3 月 21 日時点で 278 件の文書に引用されている*2ようです。これは Computer Graphics 分野では非常に多い*3引用件数です。

後続研究がたくさんありそうです。例えば同年の SGP や翌年の SIGGRAPH で発表された以下の論文は同じ研究グループからの論文で、いずれも 300 件以上の引用がある*4 *5有名研究です。

Linear Rotation-invariant Coordinates for Meshes
Yaron Lipman, Olga Sorkine, David Levin, Daniel Cohen-Or
SIGGRAPH 2005

Laplacian Surface Editing
Olga Sorkine, Daniel Cohen-Or, Yaron Lipman, Marc Alexa, Christian Roessl, Hans-Peter Seidel
SGP 2004

他にも多くの解説論文やサーベイ論文が出ていそうですが、ここでは省略します。

論文の主旨

Differential coordinates という考え方を利用して、三角形メッシュで表現された三次元形状を編集する手法を提案しています。Differential coordinates は、

In this paper we advocate the use of linear differential coordinates as means to preserve the high-frequency detail of the surface.

とあるように、サーフェスの高周波な詳細を保存するために導入されます。

さらにこの論文では編集による局所的な回転をうまく扱うための工夫をしていて、そこがこの論文の重要なところなのですが、ここでは省略します。

Differential Coordinates の概要

Differential coordinates とは座標の表現の一種で、ここではメッシュの各頂点の座標を表すために世界座標系による表現の代わりに用いることを想定しています。

日本語訳は不明です。もしかしたらまだないのかもしれません。

後述するように、Laplacian coordinates は differential coordinates の一種です。

世界座標系による表現が与えられたとき、それに対応する differential coordinates は、ある線形変換 $D$ によって得ることができます。逆に differential coordinates が与えられたとき、それに対応する世界座標系による表現は、線形変換 $D$ の逆変換 $D^{-1}$ によって得ることができます。実際には線形変換 $D$ は行列として表現されるので、逆変換を行うには線形方程式を解くことになります。

Differential Coordinates の導入

ここでは論文の Section 3: Fundamentals に当たる部分を簡単に紹介します。

まず $G=(V, E)$ を三角形メッシュとします。ここで $V$ は頂点の座標の集合、 $E$ は辺の集合です。また $j$ 番目の頂点の位置を $\mathbf{p}_j$ と表すことにします。

続いて、 $S$ を $\mathbf{p}_j \in V$ を受け取りその近似的な位置を返す写像 $S$ :

$\mathbf{p}_j \approx S(\mathbf{p}_j)$
を考えます。特にここでは、近似的な位置がその他の頂点位置の線形和で表されるものを考えます。つまり、
$\displaystyle{S(\mathbf{p}_j) = \sum_{i \in \text{supp}(j), i \neq j} \alpha_{ji} \mathbf{p}_i}$
と表される場合を考えます。ここで $\text{supp}(j)$ は $j$ 番目の頂点を近似するために使われる頂点の集合、 $\alpha_{ji}$ は $j$ 番目の頂点を近似する際に $i$ 番目の頂点位置に対して乗じられる重みを表します。 $\text{supp}(\cdot)$ と $\alpha$ を定めることで写像 $S$ は定まります。

次に、頂点位置を変換するための線形変換演算子 $D$ を

$D(\mathbf{p}_j) = \mathbf{p}_j - S(\mathbf{p}_j)$
と定義します。ベクトル $D(\mathbf{p}_j)$ を $j$ 番目の頂点の differential coordinates であると呼びます。

演算子 $D$ は線形であるため、行列形式で書くことが可能です。 $n = \left| V \right|$ とし、 $\mathbf{M}$ を $n \times n$ の行列として、

$M_{ij} = \begin{cases} 1 & (i = j) \\ - \alpha_{ij} & (j \in \text{supp}(i)) \\ 0 & (\text{otherwise}) \end{cases}$
と定義します。すると、
$\begin{bmatrix} \delta^{(x)} & \delta^{(y)} & \delta^{(z)} \end{bmatrix} = \mathbf{M} \begin{bmatrix} \mathbf{p}^{(x)} & \mathbf{p}^{(y)} & \mathbf{p}^{(z)} \end{bmatrix}$
という関係式が成立します。ただし $\delta^{(x)} \in \mathbb{R}^{n}$ は $i$ 番目の要素に $i$ 番目の頂点の differential coordinates の $x$ 座標を持つベクトルです。 $\delta^{(y)}, \delta^{(z)}$ も同様です。また、 $\mathbf{p}^{(x)} \in \mathbb{R}^{n}$ は $i$ 番目の要素に $i$ 番目の頂点の元の $x$ 座標を持つベクトルです。 $\mathbf{p}^{(y)}, \mathbf{p}^{(z)}$ も同様です。

なお、 $\left| \text{supp}(\cdot) \right|$ が小さければ $\mathbf{M}$ は疎行列になるため、元の座標から differential coordinates への変換は比較的高速です。

逆に differential coordinates の情報から元の座標情報を再構築するには、上述の線形連立方程式を解けば良い*6ことになります。その際にユーザによる形状編集の意図*7を soft constraints としてこの線形連立方程式に組み込み、それを最小二乗法で解くというのが、この手法の流れとなります。

Laplacian Coordinates との関係

Laplacian coordinates は differential coordinates の特殊形です。特に

$\displaystyle{\text{supp}(j) = \{ i \mid (j, i) \in E \}, \alpha_{ji} = \frac{1}{d_{ij}}}$
としたときの differential coordinates のことを Laplacian coordinates と呼びます。ここで $d_{ij}$ は頂点間の距離に応じた値などが使われますが、幾つかバリエーションがあります。

*1:http://dblp2.uni-trier.de/db/conf/smi/

*2:https://scholar.google.com/scholar?q=differential+coordinates+for+interactive+mesh+editing

*3:参考までに、本ブログでも散々言及している Muller らによる Position Based Dynamics という論文は現時点で 362 件でした。

*4:https://scholar.google.com/scholar?q=linear+rotation-invariant+coordinates+for+meshes

*5:https://scholar.google.com/scholar?q=laplacian+surface+editing

*6:ただし $\mathbf{M}$ が特異行列の場合は soft constraints を足さないと解が定まりません。

*7:例えば「 $i$ 番目の頂点をある位置に固定したい」などの線形な制約であれば、上述の連立方程式を拡張することで組み込むことができます。

2017-02-06

追実装: 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 の追実装を行った上で、その将来性について検討する

Projective dynamics による布のシミュレーション

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

論文読み: Projective Dynamics: Fusing Constraint Projections for Fast Simulation (SIGGRAPH 2014) - yuki-koyama's blog

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

アルゴリズム解説

準備：Implicit Euler Integration (Section 3.1)

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

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

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

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

この目的関数は、

慣性エネルギーに関する項
弾性エネルギー *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}$ に最も近い変形状態を探すための補助変数です。このように表現すると、弾性エネルギーの計算は

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

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

この処理を高速に行うために、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: Open-source C++ Library for Geometry Processing

今回実装するにあたって、クラス構造などを中心に 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 です。リアルタイム描画ではなく連番画像を出力する実装になっています。

結果

結果１：異なる解像度

この論文で 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

結果２：異なる反復回数

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

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

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

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

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

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

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

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

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

理由２：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

理由３：Projective dynamics は PBD に比べ柔軟性に欠ける

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

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

理由４：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

Sofien Bouaziz, Sebastian Martin, Tiantian Liu, Ladislav Kavan, and Mark Pauly. 2014. Projective dynamics: fusing constraint projections for fast simulation. ACM Trans. Graph. 33, 4, Article 154 (July 2014), 11 pages. DOI: http://dx.doi.org/10.1145/2601097.2601116
Isaac Chao, Ulrich Pinkall, Patrick Sanan, and Peter Schröder. 2010. A simple geometric model for elastic deformations. ACM Trans. Graph. 29, 4, Article 38 (July 2010), 6 pages. DOI: http://dx.doi.org/10.1145/1778765.1778775
Yuki Koyama, Kenshi Takayama, Nobuyuki Umetani, and Takeo Igarashi. 2012. Real-time example-based elastic deformation. In Proc. SCA '12. 19-24. DOI: http://dx.doi.org/10.2312/SCA/SCA12/019-024
Tiantian Liu, Adam W. Bargteil, James F. O'Brien, and Ladislav Kavan. 2013. Fast simulation of mass-spring systems. ACM Trans. Graph. 32, 6, Article 214 (November 2013), 7 pages. DOI: http://dx.doi.org/10.1145/2508363.2508406
Tiantian Liu, Sofien Bouaziz, Ladislav Kavan. 2016. Towards real-time simulation of hyperelastic materials. arXiv:1604.07378. URL: http://arxiv.org/abs/1604.07378
Miles Macklin and Matthias Müller. 2013. Position based fluids. ACM Trans. Graph. 32, 4, Article 104 (July 2013), 12 pages. DOI: http://dx.doi.org/10.1145/2461912.2461984
Miles Macklin, Matthias Müller, and Nuttapong Chentanez. 2016. XPBD: position-based simulation of compliant constrained dynamics. In Proc. MIG '16. 49-54. DOI: https://doi.org/10.1145/2994258.2994272
Sebastian Martin, Bernhard Thomaszewski, Eitan Grinspun, and Markus Gross. 2011. Example-based elastic materials. ACM Trans. Graph. 30, 4, Article 72 (July 2011), 8 pages. DOI: http://dx.doi.org/10.1145/2010324.1964967
Matthias Müller, Nuttapong Chentanez, Tae-Yong Kim, and Miles Macklin. 2015. Strain based dynamics. In Proc. SCA '14. 149-157. DOI: http://dx.doi.org/10.2312/sca.20141133
Matthias Müller, Bruno Heidelberger, Marcus Hennix, and John Ratcliff. 2007. Position based dynamics. J. Vis. Comun. Image Represent. 18, 2 (April 2007), 109-118. DOI: http://dx.doi.org/10.1016/j.jvcir.2007.01.005
Rahul Narain, Matthew Overby, and George E. Brown. 2016. ADMM ⊇ projective dynamics: fast simulation of general constitutive models. In Proc. SCA '16. 21-28. DOI: http://dx.doi.org/10.2312/sca.20161219
Huamin Wang, James O'Brien, and Ravi Ramamoorthi. 2010. Multi-resolution isotropic strain limiting. ACM Trans. Graph. 29, 6, Article 156 (December 2010), 10 pages. DOI: http://dx.doi.org/10.1145/1882261.1866182
Marcel Weiler, Dan Koschier, and Jan Bender. 2016. Projective fluids. In Proc. MIG '16. 79-84. DOI: https://doi.org/10.1145/2994258.2994282

*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 では制約は全てソフト制約として扱われています。

2017-01-20

最適化計算アルゴリズムCMA-ESのライブラリlibcmaesを使ってみる

ソースコード数学開発

CMA-ESについて

CMA-ES (Covariance Matrix Adaptation Evolution Strategy) は連続最適化アルゴリズムの一種です。日本語では共分散行列適応進化戦略と呼ばれます。進化戦略に基づいて、目的関数

$\displaystyle{ f : \mathbb{R}^n \rightarrow \mathbb{R} }$
を最小化する点
$\displaystyle{ \mathbf{x}^{*} = \mathop{\rm arg\,min}\limits_{\mathbf{x}}\, f(\mathbf{x}) }$
を計算するために使用されます。特徴としては

導関数を計算する必要がない (derivative-freeである)
多峰性の関数やノイズの乗った目的関数に対してもロバストである
並列性が高い

などが挙げられます。

CMA-ES - Wikipedia

このアルゴリズムは以下の論文で提案されたそうです。

N. Hansen, and A. Ostermeier.
Adapting arbitrary normal mutation distributions in evolution strategies: the covariance matrix adaptation.
Evolutionary Computation '96
https://doi.org/10.1109/ICEC.1996.542381

アルゴリズムの概要

各イテレーションで多数のサンプリング点を保持し、これらを進化計算の要領で更新していくことで最適解を発見します。進化計算の観点からは、イテレーションがgenerationに対応し、サンプリング点群がpopulationに対応しています。あるpopulationは前のgenerationのpopulationから計算される共分散行列を元に生成されます。

二次元の目的関数に対しCMA-ESを適用した様子を表した図。出典：Wikimedia Commons by Sentewolf。

libcmaesについて

CMA-ESを実装したライブラリlibcmaesがEmmanuel Benazera氏によって公開されています。ライセンスはLGPLv3です。C++11で書かれています。Python向けのバインディングもあります。

github.com

Mac OS Xでlibcmaesを動かす

下記のページに従います。

Building libcmaes on Mac OSX · beniz/libcmaes Wiki · GitHub

ます必要な依存ライブラリをbrewでインストールします。

$ brew install automake
$ brew install eigen

続いて下記コマンドを実行します。

$ git clone https://github.com/beniz/libcmaes.git
$ cd libcmaes
$ ./autogen.sh
$ ./configure --with-eigen3-include=/usr/local/include/eigen3 CXXFLAGS="-I/usr/local/include -O3" --prefix=[XXX]
$ make
$ make install

環境を変に汚したくなかったので、ここではconfigureのprefixに適当なパス (上記では [XXX] と表示) を与えました。与えたパスの配下に以下のフォルダが生成されます。

[XXX]/include
[XXX]/lib
[XXX]/bin

あとは通常のライブラリと使用方法は同じです。なおbin以下にはテストコードが含まれています。

使用例

以下のように使います。

// ヘッダファイルの読み込み
#include <libcmaes/cmaes.h>

// 目的関数の定義
libcmaes::FitFunc func = [](const double *x, const int N)
{
    double val = 0.0;
    for (int i = 0; i < N; ++ i)
    {
        val += x[i] * x[i];
    }
    return val;
};

int main(int, char**)
{
    // 探索空間の次元
    const int dimension = 10;

    // 初期解の定義
    const std::vector<double> x0(dimension, 10.0);

    // 初期の分散値
    const double sigma = 0.1;

    // 各generationにおける個体数
    const int lambda = 100;

    // 最適化計算の実行
    libcmaes::CMAParameters<> cmaparams(x0, sigma, lambda);
    libcmaes::CMASolutions cmasols = libcmaes::cmaes<>(func, cmaparams);

    // 結果の表示
    std::cout << cmasols << std::endl;

    // 返り値
    return cmasols.run_status();
}

実行結果は以下の通りです。

best solution => f-value=1.28654e-14 / fevals=13300 / sigma=0.00235696 / iter=133 / elaps=28ms / x=  7.8921e-09  1.60128e-10 -9.10801e-09 -4.36625e-08  2.52471e-08  1.12213e-09  -2.0139e-08 -1.15845e-10  2.67003e-08  9.51658e-08

ここでは10次元の関数を最適化するのに、反復回数として133回、時間として28ms必要だったことを示しています。

なぜCMA-ESか

最近読んだ論文でCMA-ESが使われているのを見て気になったので調べました。Derivative-freeなアルゴリズムは使い易いため、今後の自身の研究でも活躍する機会があるかもしれません。

Jungdam Won and Jehee Lee.
Shadow theatre: discovering human motion from a sequence of silhouettes.
SIGGRAPH 2016.
http://mrl.snu.ac.kr/research/ProjectShadowTheatre/ShadowTheatre.htm
www.youtube.com

Gaurav Bharaj, David I. W. Levin, James Tompkin, Yun Fei, Hanspeter Pfister, Wojciech Matusik, and Changxi Zheng.
Computational Design of Metallophone Contact Sounds.
SIGGRAPH Asia 2015.
http://people.seas.harvard.edu/~gaurav/papers/cdmcs_sa_2015/
www.youtube.com