中川さんのスライドを読んでいく、第三弾。

Air Meshes for Robust Collision Handling

これはデモが無いと良くわからんなぁ。 とソースコードを眺めていたら、binというフォルダに中身がある。 これ実行すればいいんじゃない？

実行出来たヽ(´ー｀)ノ

ただ早すぎて良くわからん。

AirMeshという物を生成し、そのAIr Meshが特定の制約を維持し続けるような拘束条件を用いる、という事のよう。
Air Meshというのはデモのこの黒い線で結ばれている三角形たちの事のようだ。

頂点からこれらの三角形はドロネー三角形分割というので機械的に生成出来るらしい。

これらのAir Meshの向きが変わらない、とはどういう制約になるのだろうか？

スライド51のunilateralというのは、外積の符号が変わらない、という事だろう。 だから右ネジの法則の向きが同じまま、という事で、生成された頂点が、そのなす辺をまたがない、という感じの制約になるのだろう。
体積も良くはわからんが、たぶん裏返らない、という事だろうなぁ。

これらの制約を適用すると何が起こるのだろうか？

ソースコードを見ていると、距離の拘束は別に入れているので、四角形が四角形のままなのはそちらのせいだろう。
Air Meshの制約が入るのは、メッシュの頂点が辺をまたごうとした時なので、この四角の箱がくっついて突き抜けない、というのを保証している気がする。

なるほど、衝突判定みたいなのを実装する時に、それを覆うメッシュの向きが保存される、という拘束条件でそれが表現出来る、という事か。

メッシュの品質とか最適化とか

あんまり向きが一気に変わると拘束条件を適用した状態が前の向きを再現するかは怪しいのだから、一度にboundされる拘束条件は近場では一つな事が望ましかろう。 そう考えると一気に2つ超えそうな、潰れた三角形とかがあんまり多いのは望ましくない事が考えられる。

潰れる三角形はまさに接触という状態を表す時に必要になるから全部潰れたのをなくす訳じゃないだろうが、定期的に三角形分割をやり直す方が良い、というのは納得出来る。 三角形分割は頂点同士の並びなどは変えないのだから、再分割の前に保存されていた向きは、再分割後も向きを保存し続ける限り保存されるはず。

スライド55の真ん中の絵は三角形になってないように見えるので何を意図しているのか分からない。 偶然赤と黒の点が直線になった時は直線とみなすの？それは無いような気がするが…

数式は真面目には追わないが、そんな難しい話でも無いので理解は出来る。

応用例

スライド60の応用例を見ると、布同士を重ねたいのか。 こういう時に突き抜けないように出来るってのが嬉しい訳ね。 確かにこのAir Meshという拘束条件を使えばそれは保証出来る気がする。

XPBD Position-Based Simulation of Compliant Constrained Dynamics

Air Meshgaあんまり長くなかったので同じエントリで次のも見てしまおう。

Xはextendedらしい。まぁビートXでもガンダムXでも、なんかXってつけたくなる事ってあるよな。

そしてヤング率が出て来る。懐かしいな。

スライド72〜74で式がいろいろ出てきて、アルファがゼロだとPBDと同じ、とか言うが、全然ついていけないので、まずPBDの式を見直す。

PBDの式はスライド14にあった。

これが一次のテイラー展開の式で、これが満たされるように \(\Delta p\)を決めましょう、という話だった。
sと言っていたのは、この\(\Delta p\)をどう各pに分配するか、という重みの係数みたいなものだった。

今更だが軽くここを導出しなおしてみよう。

PBDのsの導出

C(p)をコストと考えれば、gradを元にpをずらしてこれを最小化するのは、通常のgradient descentになる。 向きとしてはもっともC(p)を下げる向きはgradと同じ向きで、向きはgradの反対、となるので、 \(\Delta p\)は以下のようにおける。

ここでkは正の定数。次に、この \(\Delta p_i\)をすべてのiについて足した結果が、ちょうど一次のテイラー展開の式をゼロとするように、このkを定める。 この時点ではkはpの添字iに依存した物か。

で、ずれ具合は質量の逆数に応じて配分されるようにしたい。 そこでまず、質量の逆数をkの外に出して足し合わせる事を考える。 質量の部分で配分の割合は決定出来ているので、残りの比例定数は添字に依存しない。 これをsと置く。つまり、こうだ。

ここでsは添字によらない、正の定数。これがスライド14の(9)式だな。

あとはこのsを求める。各点に対してこの差分を足し合わせた結果が、ちょうど一次のテイラー展開の式をゼロとするようにsを決めれば良い。

これはスライド14の(8)式だな。なるほど。こうやって求まるものなのか。

点は \(p_j\) 個だけあって、それぞれの点におけるgradは違うが、いわゆる最急降下するのだから、C(p)の減らし具合に関しては絶対値だけ考えれば十分。

XPBDのラグランジュ乗数ってなんなのさ？

突然スライド72で登場するラグランジュ乗数だが、sに沿った物らしいので、sを求めてみた今ならわかるんじゃないか。

まずラグランジュ乗数と言っているのだから、何かしらのラグランジアンがある訳だよな。

何かの最小化問題を、制約条件付きで解いているはずだ。 それを考えると、制約条件はC(p)=0のような気がする。

それでは最小化したいのは？ 普通に考えたら調整にかかるエネルギーだよな。

う、これは剛性体を曲げる時のエネルギーを計算する必要がある。 引っ張るなら単純にバネと同じになるんだが、曲げるのはちょっとかったるいんだよなぁ。

ちょっとガチ変分計算になるので、ゆとりの自分には辛い。 そこでその式はとりあえず \(E(\Delta p)\)と置く事にしよう。

するとラグランジアンはたぶんこうか？

この場合、ラムダは経済学のお話で言う所のシャドウプライスに相当する訳だよな。 Cがちょっとゼロからずれた時のエネルギーコストの比例定数、になるんだっけ。 覚えてない…

気分的にはデルタpで微分してイコール0で置くと、

これにEの具体的な形を代入すると、スライド72の式(18)のようになる、という事かな。 自信は無いが。

xの更新アルゴリズム(スライド72）

合っているかはいまいち自信が無いが、あっているとして話を進めよう。

ではこのスライド72のアルゴリズム、4.1から4.4は何をやっているのか？

おそらく一回のイテレーションでだいたい必要なエネルギー最小の元でCがゼロになるようにpを動かす事は出来るのだろう。

だがPBDの時も、複数の拘束条件がある時には、イテレーションして全部の拘束条件が小さくなるような値に調整出来ると期待する、という物だった。 そこは同じだろう。

ではラグランジュ乗数を更新するとはどういう事か？

良く分からないな。 まずはラムダが特殊な値の時の振る舞いを考えてみるか。

ラムダがゼロの時は、分子はsと等しく、分母だけが少し大きくなる。 だからデルタpとしては方向は同じだが、ちょっと長さが足りない。

弾性値の分だけ少し手前で一旦止める訳だ。

ここでゼロでないラムダを足すと、どうなるのだろう？

まず、ラムダはいつも負なのかな。

では負の値のラムダをアルファと掛けて引くとデルタラムダの分子の絶対値は小さくなる…あれ？大きくなってsに近づいて行くのかと思ったが、逆なのか。

ただ、向きはgradの反対向きなので、基本的にはCの値は毎ステップで小さくなっていく。

アルファはヤング率の逆数との事なので、大きいい程ゆがみやすい訳だよな。 実際スライド70を見ても、コンクリートが一番小さく、ゴムとかが大きくなっている事から、大きい方が歪みやすい事が見て取れる。

さて、ラムダが負の値で絶対値が大きくなっていく時、何が起こるのだろうか？ 式(18)を眺めながら考えたい。

まず、ラムダの影響はアルファが大きい方が大きい。 つまり、歪みやすい物質ほど大きい。

カンニングも兼ねてデモを動かす。左右で素材が変わり、上下でイテレーションの数が…って上下はキーボードつなげてないとそんなキー無いよ！ ソフトキーである奴にしておいてくれよ… ＜ PCキーボードなしで使ってる方が珍しい

で、PBDはいまいち何に近いかわからんなぁ。数式的にはアルファの小さい、コンクリートと似た挙動になるはずなのだが、どうもそうは見えない。 むしろゴムとかに近く見える。なんでだ？

式 に戻ってラムダとアルファの事を考えよう。

式がどう出てきたかを考えると、変形にかかるエネルギーと拘束条件の兼ね合いで位置を決める、という話に見える。

ラムダとxを交互に更新していくのは、二本の微分方程式で表される関係を、ニュートンなんちゃら法とかで解いているんじゃないか？

なんだっけ。縮小写像とか使って不動点求めるんだっけ。あんま詳しくないが。

そう考えるとこれは元の二本の微分方程式を見せてくれないと解釈は難しそうだな。 そしてそうやって出る物なら、たぶんsの延長としてたどり着くのは難しそう。

ふむ。ではここまでの自分 の理解をまとめよう。

これが先のラグランジアンから出ると信じるなら、拘束条件と剛性を考慮に入れた変形を行っている事になる。 そのおかげで、何が何でも拘束条件を満たす場合に比べて不自然な動きが無くせて、剛性の効果もイテレーションを回すと厳密解に収束するんじゃないか。

それが正しいとするなら、複数の拘束条件に対して収束を待つ為には、これまで以上に多くのループを回す必要がある、という欠点がありそう。

スライドの82を見ると、計算コストは同等、とある。すると自分の理解は誤ってそうだなぁ。 各pを出さずにxが突然出せる理由があるのかなぁ。

まぁなんにせよ、このスライドで述べている事の範囲以上の事を考えてしまっているので、スライドに戻ろう。

たぶん剛性を用いて変形のエネルギーも考慮に入れて拘束条件を解く、という拡張になっていると思われる。 具体的な計算方法はスライド73の(18)式と(17)式を順番にやっていく事で計算出来る。

アルファというヤング率の逆数を入れる事が出来て、これは物理的に分かっている値を入れる事が出来るので、物理学辞典とかそういう実際の値を入れれたり、そうでないとしてもモデラーの人が入れる時に、より感覚的な値になるだろうから調整しやすかろう。

また、スライド69にはPBDではイテレーション回数やタイムステップにより剛性値が変化してしまう、という話が出ている。

PBDにどうしてそういう欠点があるかはいまいち理解出来ていないが、デモを見ていると確かにたまにぶよんぶよんなってるので、そういう欠点はあるのだろう。

XPBDでは解が変形のエネルギーなどを含んだラグランジアンから出ている（と思う）ので、イテレーションを増やしていくと厳密解、つまり与えたcompliantに従った剛性値に収束するのだと思う。 実際デモを触っていると、イテレーションの数が少ない時はぶよんぶよんしているが、増やして行くと同じような振る舞いに収束している様子が伺えるので、その改善にはなっていそう。

まとめ

中川さんのスライドを見た。 PBDがどういう物かはだいたいわかったと思う。

普通にニュートン方程式やら反発係数やらで突き抜けたりするのを考えるよりも、拘束条件をぺたぺた貼っていって制約解消系で行う、というのはエレガントに思うし、結果も安定しそうで良さそう。

さらにAir Meshのような物で複数の薄っぺらい物が折り重なるような物や、ぺったんこになるような物も扱えたり、ヤング率を各イテレーションであらわに考慮に入れるモデルで、よりリアルでありながら安定や計算コストがそれほど増えない方向に進歩しているらしい、というのは伝わってきた。

ゲームプログラミングとか全然知らない自分だが、昔物性をやっていた都合で微小体積に歪みを与えた時の全体のエネルギーとかの計算はまぁまぁ想像出来るので、素人にしては結構理解出来た方なんじゃないか。

ちゃんと知ろうと思えば個々の論文を追うべきだろうが、最近のトレンドみたいなのは伝わってきて、門外漢の自分でも「へー、そんな感じなんだ」と思える位は分かった。