物理シミュレーションの現状と未来 Physics simulation, state of the art and future ○ 長谷川晶一　東京工業大学精密工学研究所/JSTさきがけ 田崎勇一 　東京工業大学大学院情報環境学専攻 http://springhead.info/ hasevr あっと gmail.com

目次 物理シミュレーションの役割 物理シミュレーションのしくみ 計算法 安定性と誤差の解消 その他の物理シミュレーション手法 接触検出 剛体の運動．拘束力．運動方程式への拘束の組み込み 計算法 安定性と誤差の解消 前進（Explicit）積分と後退（Implicit）積分 誤差の解消法 その他の物理シミュレーション手法 Articulated Bodyのシミュレーション(Featherstone法） ペナルティ法 接触検出 大まかな判定 詳細な判定 物理シミュレータの今後 並列・専用ハード/性能と安定性 さまざまな応用

物理シミュレーションとは コンピュータの中で物理世界の物の動きを再現 物体の動きを決めるパラメータ（初期条件）を持つ 質量・摩擦係数・跳ね返り係数・位置・速度 … 物体の動きを決める法則と数値計算 運動方程式（ニュートン、オイラー、ナビエ・ストークス） 数値計算（オイラー、ルンゲクッタ、前進積分、後退積分） 未来の状態が計算で求まる 3DCGモデル：色・形 物理世界（実世界） 物理シミュレーション： 物体の運動

「ふつうの」物理シミュレーション 目的：何かを知りたい・理解したい 例： 銀河の生成 気候 観測結果にあう シミュレーション →銀河のこれま でが分かる SPH法など 気候 今日の天気の観測 →明日の天気が分かる 地球シミュレータ グリッド法 楕円銀河の光度進化[Kobayashi 2004] 大気循環モデル [地球シミュレータセンター大気・海洋シミュレーション研究グループ]

「ふつうの」物理シミュレーション 例 応力 流体 目に見えない 力を見たい 有限要素法 (FEM) 目に見えない 流れを見たい ANSYSによる構造解析の例 [ANSYS Inc] エンジン開発のための三次元流体シミュレーションの結果[NASUDA NEWS No.261]

エンタテインメントのための物理シミュレーション 目的：楽しみを感じたい エンタテインメントの対象が感じられる必要がある リアルタイム インタラクティブ 人の作用に反応 反作用を感じる そこに世界があると 感じられる リアルタイムの1/10 リアルタイム

エンタテインメントの対象と物理シミュレーション ゲーム コンピュータの中 バーチャルリアリティ コンピュータの中の世界に 入り込める 実世界指向 実世界 コンピュータが実世界を 理解・制御 コントローラ 画面 対象 計算機 計算機 没入映像 対象 インタフェース 実世界を再現 計算機 センサ 実世界のモデル 実世界を 認識・制御 ロボット 対象

コンピュータの使い方 実世界の課題をコンピュータで 旧来の使い方 実世界 の課題 新しい使い方 実世界 の課題 記録，計算，検索．．． 人間が実世界とつなぐ 新しい使い方 実世界の課題 身近な課題へ利用 身体でコンピュータを使う 　　実世界のモデル 　　人間のモデル 実世界 の課題 センサ インタフェース 計算機 実世界のモデル 人間のモデル 実世界 の課題 ロボット ディスプレイ

コンピュータの使い方 実世界志向 バーチャルリアリティ 計算機 実世界を理解・制御 するための シミュレーション コンピュータの中に 実世界を作るための シミュレーション 計算機 センサ 実世界のモデル 実世界を 認識・制御 ロボット 課題 計算機 ディスプレイ 課題 インタフェース 実世界を再現

物理シミュレーションの役割 実世界の面白さ ビデオゲーム・バーチャルリアリティ 実世界は、行動に応じて変化 多様な行動・多様な変化 ビデオゲーム・バーチャルリアリティ 行動に応じてさまざまに変化するものを作るには？ たくさんの変化の結果を用意 →　場合分けの爆発・膨大な手間と時間 物理シミュレーション →　多様でリアルな変化を自動生成 従来の新技術と同じこと 多量の2次元画像 → ３DCG ストーリの分岐 → 束構造に

バーチャルバスケット 1998

バーチャルカヌー 2005

物理シミュレーションの役割 実世界指向 実世界の便利さ、面白さをそのまま利用 計算機による実世界の認識・実世界への介入を補助 Virtual Rubik’s cube 指先の位置から手全体の位置姿勢を認識

物理シミュレーションの役割 実世界指向 実世界の認識を補助 作品Kobito 紅茶缶を画像から、こびとたちと インタラクション可能な物体として認識

まとめ 「ふつうの」物理シミュレーション →物理現象の理解予測 エンタテインメントのための物理シミュレーション →物理世界の楽しさを感じたい →物理現象の理解予測 エンタテインメントのための物理シミュレーション →物理世界の楽しさを感じたい 物理世界は楽しい 人の行動に応じて変化する 半分制御可能・半分制御不可能 物理シミュレーションは、 ゲーム・VRでは、これを再現する 実世界指向では、変化する世界を捉える

物理シミュレーションのしくみ 剛体の運動 拘束力 運動方程式への組み込み

物理シミュレーションのしくみ … m 物理法則(現実世界)は微分方程式で記述できる シミュレーションは微分方程式の数値解の計算 f x たとえば シミュレーションは微分方程式の数値解の計算 質点の運動 剛体の運動 流体の運動 運動方程式： m f 差分方程式にすると： x 順番に求めて行く： … x(0) x(Dt) x(2Dt)

剛体運動のシミュレーション 剛体 剛体のシミュレーション 硬いもの，変形しないもの 剛体でないもの 剛体だと思えるものは多い 積み木，ボール，ロボット・・・ 剛体でないもの スポンジ，粘土，水・・・ 剛体のシミュレーション 剛体だと思えるものは多い 机・椅子・箱・石… 剛体だと思えることは多い 人や動物の体の動きなど、 本当は剛体ではないけれど、剛体だと思って計算すると 運動が再現できる こんな感じに

剛体の運動 m: 質量 I: 慣性テンソル f: 外力 v: 速度 w:角速度 運動方程式 ならば，速度一定・角運動量一定

kx mg fn ft 剛体に働く力 重力→ f=mg… 定数 バネ→ f=kx… 位置に比例 拘束力 拘束力の計算が難しい 力の大きさは不明 剛体同士の位置・速度関係が決まっている 蝶番：2物体の相対位置が一定 抗力：2物体が互いに侵入しない 静止摩擦力：物体が滑らない 拘束力の計算が難しい kx mg fn ft

多体の剛体運動(Multi Body Dynamics)シミュレータ 物理シミュレータの種類 拘束力の計算法が違う 速度ベース・LCPへ定式化 もっとも良く使われている Havok Novodex ODE … 参考： Kenny Erleben: Multibody Dynamics Animation Lecture 12 Department of Computer Science Copenhagen University 多体の剛体運動(Multi Body Dynamics)シミュレータ 解析法 全自由度法 ペナルティ法 撃力法 Mirtich 96 Hahn 88 LCPによる定式化 Moore & Williams 89 加速度 ベース Baraff 89 速度ベース Stewart 96 Jean Jacques Moreau の方法 自由度削減法 Springhead1 Armstrong 79 Featherstone 83 OpenTissue ODE Springhead2

剛体の運動方程式 fc u M fe 2物体＋拘束＋外力 たくさんの剛体＋拘束＋外力 rj ri その他の外力 重力とか 速度 角速度 慣性 速度 角速度 その他の外力 重力とか 拘束力 u M fe mj Ij vj wj rj ri mi Ii vi wi

拘束の例 関節(蝶番の例、接触などは後で) w l e1 e2 e3 w1, l1 蝶番の軸e3まわりの回転以外は相対運動をしない： 連結点での相対速度・角速度 連結点に作用する拘束力・トルク w l j w2,, l2 w1, l1 w3, l3 w4, l4 w5 , l5 w6 , l6 e1 i e2 e3 拘束の座標系 速度 力 速度 力

拘束を運動方程式に組み込む w 拘束を剛体の座標系で書きたい→ w,lを書いてみる vj ,wj rj ri Jrow1 vi ,wi u 外積の行列表示 Jrow4 u w2 w3, w5 w6,も同様に書ける

拘束を運動方程式に組み込む vj ,wj rj ri vi ,wi fj ,tj rj ri fi ,ti w l u 剛体座標の速度・角速度 l fj ,tj rj ri fi ,ti fc 剛体座標での 力・トルク l 拘束座標での力・トルク

拘束の組み込み 拘束を考えやすいように、運動方程式を変形 拘束条件と連立させる(蝶番の例) 拘束： 運動方程式に代入： 連立方程式を解いて求める

シミュレーションの手順 拘束力l を計算 剛体の速度を更新 剛体の位置を更新 計算法のところで説明 全剛体の 位置姿勢 位置 姿勢 Quaternion 速度を Quaternionの微分に 変換する行列

接触の拘束 接触 お互いに侵入しない ＝力が働いて止まる 又は 力は働かず離れる： 静止摩擦か動摩擦： 速度 トルクは0： 力 速度 力 接触点での相対速度 接触点に作用する接触力 速度 力 速度 力 速度 力

拘束の組み込み(接触の例） 拘束を考えやすいように変形した運動方程式 拘束条件と連立させる 速度 力 速度 力 一見、 1本の式に変数が2つあるように 見えるけれど、 実はどちらか片方しか動かないので解ける

まとめ 速度ベースの拘束法・LCPへ定式化する方法 が良く使われている LCPに定式化して拘束力を求め、積分を繰り返す 拘束 剛体の運動方程式をヤコビアンJで座標変換 → wとlの運動方程式 拘束を相対速度wと拘束力lについての式に定式化 wとlの運動方程式に代入→LCPになる LCPを解いてlを求め、剛体の速度・位置を更新 拘束 関節は直線、接触は折れ線 ２変数あるように見えて実は１変数

連立方程式の解法 LCPの近似解法 誤差の解消法 計算法 連立方程式の解法 LCPの近似解法 誤差の解消法

連立方程式の解き方 拘束を組み込んだ運動方程式 LCPを解く 線形相補問題 (LCP:Linear Complementary Problem) LCPを解く ピボッティング法　・・・　厳密．低速 反復解法　・・・　高速．精度は反復回数に依る 行列分割法 ヤコビ法 ガウス・ザイデル法,SOR法 ニュートン法 式1本で変数が2つ

ガウスザイデル法・SOR法 行列分割法 巨大な連立方程式の近似解を少ない計算で解く 方法のひとつ Dは対角行列のような簡単な行列にする。 DをAの対角成分とし、残りをFとしたのがヤコビ法 それを工夫したのがガウスザイデル法 もし、漸化式 が収束するなら なので、求めたい となる。 Aが正定値対称行列なら収束する 細長いもの・慣性テンソルが小さいものは 収束しにくい

ガウスザイデル法・SOR法 ガウスザイデル法 SOR法 ガウスザイデル法をちょっとだけ加速 早く収束することもある ヤコビ法：lを一度に更新 更新済み 更新中 未更新 ガウスザイデル法 1.6倍加速

l1を更新するたびに、 l1>0をチェック 負になったら l1=0に設定 LCPとガウスザイデル法 ガウスザイデル法は1行ずつ更新 接触力をうまく計算できる 更新済み 更新中 未更新 　お互いに侵入しない＝力が働いて止まる 又は 力は働かず離れる： 　静止摩擦か動摩擦： l1を更新するたびに、 l1>0をチェック 負になったら l1=0に設定 更新したl1を使って -ml1 <l2< ml1をチェック はみ出したら、l1=+-m l1でクリップ

シミュレーション例

解析法とペナルティ法の計算の比較 どちらも繰り返し計算で接触力を求めるが： 解析法は，最初から接触力同士の相互作用を考慮 少ない回数で解にたどり着く可能性がある 解析法（ガウスザイデル法） ペナルティ法

前進（Explicit）積分と後退（Implicit）積分 バネダンパの例 拘束を満たす時刻と安定性 誤差の解消法 安定性と誤差の解消法 前進（Explicit）積分と後退（Implicit）積分 バネダンパの例 拘束を満たす時刻と安定性 誤差の解消法 この項については，以下で発表をしています： 田崎勇一, 長谷川晶一, “拘束法の動力学シミュレータのための安定なバネダンパモデル”, 情報処理学会 グラフィクスとCAD研究会第124回研究発表会資料, 2006.8

K B 安定性と積分の方法 前進（Explicit）積分と後退（Implicit）積分 バネ・ダンパの運動方程式 m Explicit Euler 法 バネ・ダンパ係数やステップ幅を大きくすると不安定となる

バネ・ダンパの計算法 （Explicit Euler法） なぜ不安定か？ 安定となるための バネ・ダンパ係数の領域 現在時刻の位置・速度から現在時刻の力を計算 →　積分誤差による不安定化

バネ・ダンパの計算法 （Implicit Euler法） 次の時刻の位置・速度から現在時刻の力を計算 全域で安定

拘束を満たす時刻と安定性 加速度ベース [baraff 89など] 積分誤差の累積によるドリフト 速度更新 位置更新 運動方程式 解析法 全自由度法 LCP 加速度 ベース Baraff 89 速度ベース Stewart 96 Jean Jacques Moreau の方法 自由度削減法 Armstrong 79 Featherstone 83 OpenTissue ODE Springhead2 加速度ベース　[baraff 89など] 運動方程式 拘束式 （加速度レベル） この瞬間[t]に、拘束を満たす 連立 速度更新 位置更新 積分誤差の累積によるドリフト

拘束運動の計算法 速度ベース [Stewart96など] 積分（速度更新）を考慮して拘束力を計算→ドリフトが少ない 運動方程式 連立 解析法 全自由度法 LCP 加速度 ベース Baraff 89 速度ベース Stewart 96 Jean Jacques Moreau の方法 自由度削減法 Armstrong 79 Featherstone 83 OpenTissue ODE Springhead2 速度ベース [Stewart96など] 運動方程式 連立 速度更新 拘束式 （速度レベル） 次の時刻[t+1]に拘束を満たす 位置更新 積分（速度更新）を考慮して拘束力を計算→ドリフトが少ない

拘束は、w[t+1]についての拘束になっている 拘束の組み込み 拘束を考えやすいように、運動方程式を変形 拘束条件と連立させる(蝶番の例) ここが、w[t+1]なので、 拘束： 拘束は、w[t+1]についての拘束になっている 運動方程式に代入： 連立方程式を解いて求める

Implicitのバネダンパ 時刻[t+1]にバネダンパの制約を満たす LCPに組み込むことができる 関節にバネダンパー 並進のバネダンパー 関節角速度 関節トルク バネダンパ力 伸びの速度 関節角度 伸び（定数）

安定性の評価 剛体5つ 蝶番5つ，4つにバネダンパ バネダンパ1つ 重力 このバネ・ダンパについて Implicit版とExplicit版を比較

誤差の解消法 誤差がたまる原因は…もう一度拘束を考えてみる 誤差の解消法 vj ,wj rj ri 速度を変化させる方法（ペナルティ法） ODE, 誤差が残る(隙間が開いたりする)が自然で安定 位置を変化させる方法（位置についての解析法） OpenTissue, 誤差が見えないが、不安定になることがある vj ,wj rj ri vi ,wi 関節位置での相対速度がwi[t+1]=0 となるように、拘束力lを計算している 速度を積分した相対位置には、誤差がたまる →誤差を解消する必要あり p,w

誤差が減るような速度になるような拘束力l’が求まる 誤差の解消法１：速度を変化させる方法 e 拘束＋運動方程式 に速度の修正を加える 誤差が減るような速度になるような拘束力l’が求まる

誤差の解消法２：位置についての解析法 運動方程式から作った速度の更新式 の位置版を考える。 e 位置版に、拘束条件の変形をする 速度版 ガウスザイデル法でlを求め、s’を求める

誤差の解消法：どちらがよいか 速度を変化 位置の解析法 誤差はすぐには解消しない ペナルティ法的 わずかな式変形 計算はとても少ない 誤差はすぐには解消しない 　ペナルティ法的 わずかな式変形 　計算はとても少ない 速度が変化する 位置の解析法 毎回、誤差を0にしてくれる 　解析法的 Al+bのAは再利用できるが 　あとは計算が必要 　ガウスザイデルをもう1つ解く 速度に影響がない 誤差 誤差 ステップ ステップ

誤差の解消法：どちらがよいか 強い外力の影響 強い外力feに対抗する 関節の拘束力l ガウス・ザイデル法の 打ち切りのため、 | l | < | fe | となり、剛体はだんだん加速。 速度変化で解消の場合 位置の解析法で解消の場合 e fe l w で位置sを補正するが、速度wは補正なし→加速し続ける 速度変化は必須、位置の補正は誤差を隠してくれる

ガウスザイデル法の誤差 計算のイメージ どちらから近づくか 1.0 1.2 lprev 0.8 内側（原点側）から近づく →真の値より小さな値 →安定、ダンパ気味 外側から近づく →真の値より大きな値 →不安定 １つ前のシミュレーションで求めた lprevからスタートすると収束が早い が、外側から近づく危険あり 　→　不安定になりやすい（例） 1.0 1.2 lprev 0.8

その他の物理シミュレーション手法 フェザーストーン法 ペナルティ法

Articulated Bodyのシミュレーション 解析法 全自由度法 LCP 加速度 ベース Baraff 89 速度ベース Stewart 96 Jean Jacques Moreau の方法 自由度削減法 Armstrong 79 Featherstone 83 OpenTissue ODE Springhead2 Articulated Body 多数の剛体を関節でつないだもの 動物や人間の体，ロボットなど 普通にLCPでも解けるが… 自由度削減法 可動関節だけをシミュレーション 剛体が輪になっていない＝木構造 で高速な手法 輪になっている例

解析法で求めると．．． l4 l3 l2 l5 l1 l6 からlを解いて、剛体の速度uを更新 uから剛体の姿勢sを更新

Featherstoneの方法 状態変数 リンクの位置・速度など他の値は持たない 根のリンクの速度・角速度，位置・姿勢 各関節の可動部分の角度，角速度 リンクの位置・速度など他の値は持たない 関節が外れることはない リンクi 根のリンク リンクi と，より葉側のリンク 根側のリンク 葉のリンク

関節の拘束力の計算 Featherstoneの方法 木構造のリンク用 リンクiの運動方程式は， のように，書ける． 　のように，書ける． Ziはリンクiと葉側のリンクの 姿勢・速度で決まる Ii , Zi を葉から順に計算 根側から，v, wを計算 Ii: リンクIと葉側のリンクの慣性項 Zi: 外力，重力，コリオリ力による項 Ii,Ziは関節加速度によらない リンクi 根のリンク リンクi と，より葉側のリンク 根側のリンク 葉のリンク

FeatherstoneのLCPへの組み込み 解析法 全自由度法 LCP 加速度 ベース Baraff 89 速度ベース Stewart 96 Jean Jacques Moreau の方法 自由度削減法 Armstrong 79 Featherstone 83 OpenTissue ODE Springhead2 LCPによる定式化では、 剛体の質量・慣性を並べて運動方程式を作った これにArticulated Bodyを加えると IA wjoint wjoint wjoint wjoint vbase wbase wjoint wjoint 参考： S. Redon et al. "Adaptive Dynamics of Articulated Bodies" SIGGRAPH 2005

解析法 ここまでで取り上げた方式 利点 動作が安定 低速更新が可能 欠点 1ステップの計算量はある程度多い. ⇒低速更新 安定だが精度が悪い 全自由度法 ペナルティ法 撃力法 LCPによる定式化 速度ベース Jean Jacques Moreau の方法 自由度 削減法 ここまでで取り上げた方式 利点 動作が安定 低速更新が可能 欠点 1ステップの計算量はある程度多い. ⇒低速更新 安定だが精度が悪い David baraff introduce analytical methods to dynamics simulators for graphics since 1989. Analytical method solves equation of motions and constraints as a linear equation system. So the number of the equations are depends on the number of objects and contacts. It takes computation order of third power of n to solve the equation system.

. . ペナルティ法 拘束を満たす向きに適当な力を加える 侵入量 d，侵入量の速度 d 滑り量 l，滑り速度l バネ ダンパ バネ ダンパ 撃力法 解析法 全自由度法 LCPによる定式化 速度ベース Jean Jacques Moreau の方法 自由度 削減法 ペナルティ法 拘束を満たす向きに適当な力を加える 侵入量 d，侵入量の速度 d . バネ ダンパ . 滑り量 l，滑り速度l バネ ダンパ 利点 1ステップの計算が速い. O(n) ⇒高速更新可能 高速更新⇒高精度 Featherstone法と簡単に組み合わせられる 欠点 Dｔを小さくしなければならない 安定性が低い Before analytical methods, penalty methods has been used. Penalty method solves constraints iteratively using spring-damper model. The constraints force is defined by the length of spring and relative velocity of damper. In this method, forces are directly defined. So the computation order is just n.

接触体積の形状を考えたペナルティ法 広い接触面があるとき バネダンパモデルはどこに置けばよいか？ ? この絵が見づらい 最侵入点にした場合

? 広い接触面の扱い バネダンパモデルはどこに置けばよいか？ 頂点にした場合 Top view Putting on each vertex can be a candidate. But it does not simulate friction force well. This is a top view of the scene. When the object rotates, a friction force is generated. The friction force should caused by whole contact area. So the friction force from these 4 vertices are not correct.

? 広い接触面の扱い バネダンパモデルはどこに置けばよいか？ 多点でサンプリング うまく働くが， 計算時間がかかりすぎる Putting many spring-damper model can be a solution. But, it takes much computation time and memories. 多点でサンプリング

! 広い接触面の扱い 分散バネ・ダンパモデル 3角形ごとに分散バネ・ダンパモデルからの力を積分 So, rather than putting many models, we have better to suppose distributed spring-damper model on the contact area. And integrate forces from the model for each triangles. So, we choose distributed models. And it realizes fast and accurate calculation of contact forces. 分散バネ・ダンパモデル

処理の手順 接触力の計算 接触体積内の1点と法線を求める（GJK) 接触体積の形状を求める(Muller and Preparata) 接触面にかかった力を積分する Next, I’ll explain the method to calculate the contact force. To find the contact force, our simulator does this: 1st: Find contact point and normal. 2nd: Find the shape of the contact volume. 3rd: Integrate penalty force over the contact area. I’ll explain each.

接触解析 接触部分= 2つの凸多面体の交差部分. D. E. Muller and F.P.Preparata: “Finding the intersection of two convex” (1978) 2つの凸多面体と交差部分上の1点が与えられたとき， 交差部分の形状を求める Because objects are represented by convex polyhedrons, the contact part of two objects is intersection of two convex polyhedrons. Muller and Preparata propose a paper entitled “Finding the intersection of two convex” in 1978. It finds intersection of two convex polyhedrons from the faces of the two convex polyhedrons and a point in the intersection. I’ll explain the algorithm briefly.

分散バネ・ダンパモデルからの力を3角形ごとに積分 力の積分 抗力 動摩擦力 最大静止摩擦力 Penalty forces and dynamic friction force can be integrate for each triangle. However, for static friction force, we need distributed spring-damper model for sliding. The slide motion of two solid can be represented by translation and rotation. So, the distributed springs can be represented by two models. So, we need just two spring-damper models for static friction. It decreases computation time drastically. 分散バネ・ダンパモデルからの力を3角形ごとに積分

3角形ごとの積分 p3 p1 p2 h3 h2 p3 抗力バネモデルによる力 抗力バネモデルによるトルク

フォースフィードバックへの応用例 高速更新なので、フォースフィードバックの感覚が良い

ペナルティ法と解析法の拘束力と動作の比較 次のステップで拘束を満たすlを 繰り返し計算で求める ペナルティ法 拘束違反を減らす向きのlを 毎ステップ加える 位置の 誤差 l20 lp 速度の 誤差 l\

大まかな判定 バウンディングボックス ソート 詳細な判定 GJK 接触形状 接触検出 大まかな判定 バウンディングボックス ソート 詳細な判定 GJK 接触形状

接触検出 接触の拘束を考えるためには、接触検出が必要 高速化のための工夫 まず大まかな判定をして、明らかに接触してないものを候補からはずす。 その後詳細な判定をする。

接触検出と計算速度 接触力を求める→高速な衝突判定が必要 階層化 境界（Bounding）の例 K-DOP 球 AABB ぴったりフィット 簡単な境界で囲む 階層的に境界を作る 大まかに判定してから，詳細な判定をする 境界（Bounding）の例 K-DOP 球 AABB ぴったりフィット 判定簡単

階層化と計算時間 log2n 階層化すると　　　　　　回の検出で終わる log2n n

アルゴリズムの速度 O(n) O(log n) O(1) a を検索 ランダム ソートされている場合 頭から見ていく きちんと並んでいる場合 A d h j J L f C g M a l m D b p E H U V s t w Q R W Y x z O Z O(n) 頭から見ていく A C D E H J L M O Q R U V W Y Z a b d f g h j l m p s t w x z あたりをつけてみていく O(log n) A B C D E F G H I J K L M N O P Q R S T U V W X Y Z a b c d e どこにあるか分かる O(1)

アルゴリズムの速度 対象に対する知識があるほど速く出来る ランダム ソートされている場合 きちんと並んでいる場合 知識なし a = = ? しか分からない場合 ソートされている場合 左＜右 ＜演算が定義されている きちんと並んでいる場合 右－左＝3 何がどこにあるか分かっている

アルゴリズムの速度（ソート） 総当りソート： O(n2) クイックソート： O(n log n) バケツソート： O(n) 1 2 5 7 6 9 3 5 7 3 6 9 1 2 <5を前へ n回比較 3 1 2 5 7 6 9 log n回 3 1 2 <3を前へ 7 6 9 <7を前へ n回比較 5 7 3 6 9 1 2 1 2 3 5 6 7 9 1 2 3 4 5 6 7 8 9 n回移動 1回

O(n2 ) 運動する物の接触検出 物体が動き回る場合，毎回境界の更新が必要 １ ２ ４ ３ ５ ６ n n 境界の更新に時間がかかる 総当りも時間がかかる 1-2 1-3 2-3 1-4 2-4 3-4 1-5 2-5 3-5 4-5 1-6 2-6 3-6 4-6 5-6 １ ２ ４ ３ ５ ６ n n O(n2 )

O(n log n) 運動する物の接触検出 ソート＝境界を作りながら判定 １ ２ ５ ４ ３ ６ １ ２ ４ ４ ５ ３ ６ 1より左？ 1-2 1-3 1-4 1-5 1-6 ５ ４ ３ ６ 1より右？ 1-3 1-4 1-6 2より左？ 2-5 2-4 １ ２ ４ ４ ５ ３ ６ 4より左？ 2-4 要判定 1-4 の判定 O(n log n)

接触検出 中身の判定 凸形状 非凸形状 中身は 三角形？ 多面体？ 凸多面体 凸形状 距離が極小となる点が1点 中身は 三角形？ 多面体？ 凸多面体 凸形状 距離が極小となる点が1点 最近傍点が簡単に求まる GJK algorithm E. G. Gilbert, D. W. Johnson and S. S. Keerthi A Fast Procedure for Computing the Distance between Complex Objects in Three-Dimensional Space (1988) ほかにLin-Canny algorithm というのもあります Ming Lin, John Canny A fast algorithm for incremental distance calculation (1991) 凸形状 非凸形状

GJK 凸形状A上の点から，凸形状B上の点へのベクトルを 原点を始点に並べる（Minkowski sum）と ベクトルの終点の集合も凸形状になる 元の図形 Minkowski Sum をとったもの

GJK ３ ４ １ ２ OVi Wi V0 : 凸形状内の任意の点 Wi :OViとOWiの内積が最小の点(support point) Vi :三角形Wi-2 Wi-1 Wi内の点で原点に一番近い点 Vi = s Wi-2 +t Wi-1 +(1-s-t) Wi W1 W2 W3 V 元の図形上で、最近傍点が求まる

速度を考慮したGJK 衝突の瞬間（Time of Impact TOI）と 衝突の位置・法線が分かる． 接触点 もとの図形 Minkowski Sum

GJKを使った接触点の計算 １ ３ + + + ２ ４ + + 参照：　Continuous Collision Detection of General Convex Objects under Translation [G. van den Bergen]

GJKの結果の解釈 法線 そのまま実世界の法線 接触点 Support point を内分した点 法線 そのまま実世界の法線 接触点 Support point を内分した点 →実世界の Support points を内分して求める． これを求める→実世界の頂点の引き算 どれを使ったか覚えておけば元に戻せる

GJKの結果の解釈２ 重なっていると，戻せない 元の形状でのの3つの点は，変換した形状では1点 →どこだか分からない． 元の形状 Minkowski Sum 元の形状でのの3つの点は，変換した形状では1点 →どこだか分からない．

接触解析 衝突部分の形状（切り口）を求める方法 切り口の頂点に接触の拘束を考える D. E. Muller and F.P.Preparata: “Finding the intersection of two convex” (1978) 凸形状２つの交差部分の形を求める 共有点１点が必要（GJKで求められる） 切り口の形、頂点を知ることができる Because objects are represented by convex polyhedrons, the contact part of two objects is intersection of two convex polyhedrons. Muller and Preparata propose a paper entitled “Finding the intersection of two convex” in 1978. It finds intersection of two convex polyhedrons from the faces of the two convex polyhedrons and a point in the intersection. I’ll explain the algorithm briefly. 切り口 上から見た図 断面図

接触解析(2) ２つの凸形状の交差部 半空間表現 Muller’s algorithm represents two convex by half space. Half space representation is common area of half spaces, which represented by planes. So, the intersection set of the two half space representations results in the intersection of the two convex polyhedrons. But, the representation has redundant planes. So we have to find the minimum set of the planes which represents the intersection. 半空間表現

接触解析(3) ２つの凸形状の交差部(2) n log n 半空間表現 双対変換 交差部の頂点 双対変換 Quick hull To find the minimum set, the algorithm uses dual transformation. Dual transformation transform a vertex into a plane and a plane into a vertex. The dual transformation of a dual transformation is identical transformation. The problem of the minimum set of the planes results in a convex hull problem, which find the minimum convex from the vertices. There is fast algorithm to solve convex hull problem. We use quick hull algorithm for this. After finding the convex hull, Muller’s algorithm uses dual transform again. Then we get the vertices and faces of the intersection. Quick hull n log n

並列化・ハードウェア化 スピードと安定性 いろいろな利用 物理シミュレーションのこれから 並列化・ハードウェア化 スピードと安定性 いろいろな利用

並列化について 物理シミュレーションは たくさんのものに同じ法則を適用 →並列化に向いている 専用計算機向きか？ 粗い判定： n2の並列計算も可能 GJK・接触解析： 接触ペアの数だけ並列計算 ガウスザイデル： 行列の1行ごとに並列計算 専用計算機向きか？ GRAPE 天体シミュレーション用 (1991　15GFPS ＠ i486の時代) 重力計算が主 → 専用ハードが作りやすい 剛体・流体…は？ データが複雑（特に接触検出・解析） 対象がさまざま →専用チップ化にはあまり向かない CELL/PhysX は並列汎用計算機。 GPGPUはソートなどあと一歩。

スピードと安定性 剛体はスピードより安定性 そろそろ数は十分 どんなシミュレータも発振する 制御の安定性の保証がほしい PD制御発振の例 Wrotek et al. 2006 Havok FXの例

スピードと安定性 流体はまだまだスピードがたりない 2次元の波面 → 十分リアルタイム 3次元のメッシュ → かなり厳しい 2次元の波面 → 十分リアルタイム 3次元のメッシュ → かなり厳しい 3次元の粒子 → たくさんの量は扱えない 組み合わせが必要かもしれない A Fluid Resistance Map Method for Real-time Haptic Interaction with Fluids [Y. Dobashi et al.　2006] Efficient Simulation of Large Bodies of Water by Coupling Two- and Three-Dimensional Techniques [G. Irving et al. 2006] R 力覚フィードバック 圧力場・流速 (FRM :Fluid Resistance Map) 3次元流体シミュレーション 水の速度 リアルタイム 事前計算 パドルの速度 リアルタイム波シミュレーション

物理シミュレーションの利用 リアルな動きの生成に利用 実世界とバーチャル世界をつなぐ ゲームの物理 キャラクタモーション 身体のシミュレーション 実世界とバーチャル世界をつなぐ 実世界の認識の補助 実世界で再現可能なバーチャル世界を作る

キャラクタモーション Oshita 2001 Liu et al.2005 軌道全体を最適化 　　軌道全体を最適化 　○ 高度な最適化 　 　△ リアルタイム生成は困難 Zordan et al. 2005,　 Zordan and Hodgins 2002, Natural Motion Endorphin コントローラ＋シミュレーション．　関節角度をPD制御 　○ リアルタイム生成 　△ コントローラの作成 PD制御で加速度を決める ← トルク/ZMPの制約 　○ PDの調整が簡単。 　△外力が加わった場合の動作 トルクではなく、加速度

キャラクタモーション モーションキャプチャデータの変換 関節トルク・床反力の計算 人体モデルのマーカ位置を キャプチャした データ位置に PD制御 関節角を記録 関節トルク・床反力の計算 関節をPD制御して動作を再生 関節と床から働く力を計算

キャラクタモーション Hasegawa et.al 2005 関節と手先位置をPD制御 物理シミュレータ 力覚インタフェース The World 力フィード バック 手の位置 物理シミュレータ Virtual human’s mind 行動ルール Motion Controller Reaching motion Default posture 予測用シミュレータ State transition Wait Guard Attack 力覚インタフェース

身体のシミュレーション 筋骨格系＋神経系のシミュレーション M. Otake and Y. Nakamura 2005

ロボット・作業のシミュレーション レスキュー ロボカップ 機構シミュレーション

まとめ エンタテインメントのための物理シミュレーション 物理シミュレーションのしくみ 物理シミュレータの今後 人が物理を感じられるシミュレーション リアルタイム・インタラクティブ 物理シミュレーションのしくみ 定式化、LCPの近似解法 安定性と誤差の解消 Featherstone法・ペナルティ法 接触検出（GJK, 接触部分の形状） 物理シミュレータの今後 並列・専用ハード/スピードと安定性 さまざまな応用